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ABSTRACT 

' The peculiar velocity of the intergalactic gas responsible for the cosmic 21cm back- 

I ground from the epoch of reionization (EOR) and beyond introduces an anisotropy in the 

>«*^ , three-dimensional power spectrum of brightness temperature fluctuations. Measurement of 

(~| ■ this anisotropy by future 21cm surveys is a promising tool for separating cosmology from 

21cm astrophysics. However, previous attempts to model the signal have often neglected pe- 
culiar velocity or only approximated it crudely. This paper re-examines the effects of peculiar 
velocity on the 21cm signal in detail, improving upon past treatment and addressing several 
issues for the first time. (1) We show that even the angle-averaged ^ov^tr: spectrum, P{k), 
^ ■ is affected significantly by the peculiar velocity. (2) We re-derive the brightness temperature 

dependence on atomic hydrogen density, spin temperature, peculiar velocity and its gradient, 
and redshift, to clarify the roles of thermal vs. velocity broadening and finite optical depth. 
(3) We show that properly accounting for finite optical depth eliminates the unphysical diver- 
gence of the 21cm brightness temperature in overdense regions of the IGM found by previous 
work that employed the usual optically-thin approximation. (4) We find that the approxima- 
' tion made previously to circumvent the diverging brightness temperature problem by capping 

\ the velocity-gradient can misestimate the power spectrum on all scales. (5) We further show 

that the observed power spectrum in redshift-space remains finite even in the optically-thin 
approximation if one properly accounts for the redshift-space distortion. However, results that 
take full account of finite optical depth show that this approximation is only accurate in the 
limit of high spin temperature. (6) We also show that the linear theory for redshift-space dis- 
tortion widely employed to predict the 21cm power spectrum results in a ~ 30% error in 
the observationally relevant wavenumber range fc ^ 0.1 — 1 /i/Mpc, when strong ionization 
fluctuations exist (e.g. at the 50% ionized epoch). We derive an alternative, quasi-linear for- 
\ mulation which improves upon the accuracy of the linear theory. (7) We describe and test two 

. numerical schemes to calculate the 21cm signal from reionization simulations to incorporate 

peculiar velocity effects in the optically-thin approximation accurately, by real- to redshift- 
space re-mapping of the H I density. One is particle-based, the other grid-based, and while the 
former is most accurate, we demonstrate that the latter is computationally more efficient and 
can be optimized so as to achieve sufficient accuracy. 

Key words: Cosmology: theory-reionization- physical data and processes: radiative 
transfer- methods: analytical-numerical- galaxies: intergalactic medium 
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1 INTRODUCTION 

Neutral hydrogen atoms in the intergalactic medium (IGM) at high 
redshift produce a diffuse background of redshifted 21cm radia- 
* Email: ymao@astro.as.utexas.edu tiori which encodes information about the physical conditions in the 

t Email: shapiro@astro.as.utexas.edu early universe during and before the epoch of reionization (EOR, 
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z > 6). Three-dimensional mapping of this 21cm background 
(a.k.a. 21cm tomography) has recently been proposed as a promis- 
ing cosmological probe. In principle, it has greater potential than 
the cosmic microwave background (CMB) since it can map most 
of our horizon volume, thus pro viding unprecedented cosmologi- 
cal information jMao et al.ll2008h . 

The next few decades promise to become a golden age for 
21cm tomography, with about a half-dozen experiments already 
proposed or underway for measuring the 21cm background from 
the EOR, including the upcoming first generation such as 21Cm4E 
MW^E LOFAkETgMRtQ, and PAPEfQ, an d the next ge neration 
such as SK^ and the Omn; j-copePll lTegmark & Zaldarriagal2009i 
These telescopes will measure the 21cm signal either statis- 
tically (first generation telescopes) or by precise imaging and map 
making (next generation telescopes). 

Observations will measure the power spectra of 21cm bright- 
ness temperature fluctuations from the EOR. The information that 
21cm power spectra encode is twofold. First, cosmic reionization 
leaves its imprint, such as the size distribution of the H II region, on 
21cm power spectra. Since the topology and geometry of ionized 
bub bles is sensitive to th e properties of the ionizing sources (see, 
e.g.. lFriedrich et al.ll20r l|), we can learn about the ionizing sources 
from 21cm power spectra. For example, we can distinguish models 
with only high-mass atomic cooling sources from models with both 
high-mass and s elf -regulated low-mass atomic cooling sources 
jlliev et alJl20Tlh . Second, 21cm power spectra are also sensitive 
to cosmological parameters because the latter determine the mat- 
ter density fluctuations at high redshifts. The precision with which 
21cm tomography can constrain cosmological parameters has been 
forecast in several studies. Some of these consider mapping dif fuse 
hydrogen in the IGM before and during t he EOR CMcOuinn et al.l 
20061: iBowman. Morales & Hewittll2007l: ISantos & Coorav|l200a 



Mao et al.1 l2008l : iBarger et al.l l2009l : lAdshead et al.l 1201 1[) , oth- 
ers mapping neutral hydrogen in galactic halos after reio nization 
l lWvithe. Loeb & Geil|[2008l ; IVisbal. Loeb & Wvithell2009l) . These 
studies show that cosmological constraints based on CMB mea- 
surements can be significantly improved if combined with 21cm 
measurements. In addition, it has been demonstrated in the litera- 
ture that 21cm power spectra can also constrain many cosmolog- 
ical models beyond the vanilla ACDM model, e.g., spatial cur- 
vature and the ru nning of the spectra of primordial s calar den- 
sity pert urbations (iMao et al.l 20081: Barger et al. |l2009l) . neutrino 
masses iMao et alj boOSi : iPrit chard & Pierpaoli '2008"), compen- 
sated isocurvature perturbations (Gordon & Pritchard 2009) , pri - 
mordial non-Gaussia n density perturbations djoudaki et al.ll201lh . 
cosmic string wakes jBrandenberger et al.l^20 10l), and anisotropic 
matter density fluctuations (Hernandez & Holderl l201 ll) . 

In view of this promise which observations of 21cm power 
spectra hold for testing and constraining cosmological and astro- 
physical models, further progress is required to ensure that pre- 
dictions are accurate enough to fulfill this promise. This accuracy 
depends not only on the realistic astrophysical modeling of reion- 
ization and the H I spin temperature, but also on the methods used 



http : / /21 cma -bao ■ ac ■ cn^ 



http : / /ww w .haystack . m it ■ edu/ast/arrays/mwa/| 
http : / / www . lof ar . orgj 
http : / /gmrt . ncra . tif r . res . in] 
http : / /astro .berkeley ■ edu/ -dbacker/eorT] 



to extract the 21cm signal from simulations of the EOR. We fo- 
cus here on this 21cm methodology issue, and leave aside the issue 
of the accuracy of the underlying reionization models and simula- 
tions. For this purpose, we will make use of the results of a recent 
reionization simulation of our own, based upon a radiative transfer 
calculation combined with a high-resolution N-body simulation of 
ACDM. While this simulation represents the current state-of-the- 
art in large-scale reionization simulations, it will serve here only 
as our illustrative testbed. The accuracy and realism of the simu- 
lation, itself, is not our concern here, as we focus, instead, on the 
accuracy of our method for extracting the 21cm signal from such 
simulations. 

All observations will give the redshifted 21cm signal in ob- 
server redshift-space, where the frequency not only depends on 
the cosmological redshift, but also on the peculiar velocity of the 
IGM. Howe ver, most theoretical endeavors, both i n anal ytical mod- 
ellin g (e.g., iFurlanetto. Zaldarriaga & Hemguist 2004: Ili ey et al] 
12002^, semi-numerical (e.g., Alvarez et al..,2009i: .Zah n et al.'2007, 



201 1) and in numerica l simulations (e.g..lShapiro et al. 2006. 200| , 
iMellema et al. 2006b|, llliev et alj 12008 d : iMcOuinn et alj |2007 ; 



iTrac & Cen 2007), have focused on predicting the statistics of 
the 21cm signal (e.g. the power spectrum of brightness tempera- 
ture fluctuations) without taking peculiar velocities into account. 
On the other hand, peculiar velocities will influence the 21cm 
brightness temperature significantly as was for example shown by 
iMellema et a l. (2006b). 

In the linear regime, the effects of peculiar velocit i es hav e 



been studied an alytically by 'Bharadwai, Nath & Sethi' ll200Ih. 



iBharadwai &Mi (2004), Barkana & Loeb (2005) and Wan g & Hul 
( l2006h . It has been shown that, in this regime, it is possible to 
separate the contributions to the brightness temperature fluctuation 
statistics from the patchiness of r eionization and the cosm ological 
density fluctuations, respectively dSarkana & LoeblboOSh . This, it 
is hoped, would make it possible to use 2Icm measurements to 
solve for cosmological parameters. However, the effects of non- 
linearity remain largely unexplorecfl There are two kinds of non- 
linearity that may contribute, one associated with the gravitational 
growth of matter density and velocity perturbations, the other due 
to ionization patchiness. The line ar th eory formula for the 21cm 
redshift-space power spectrum ( Barkana & Loebl2005i) widely em- 
ployed in the literature was derived under the assumption that not 
only the matter density and velocity fluctuations are linear, but so 
are the ionization fluctuations. The latter assumption clearly breaks 
down on the scale of the size of the H II region. We shall inves- 
tigate here the accuracy of this linear theory formula, particularly 
for the wavenumbers that are expected to be probed by current and 
future 21cm surveys of the EOR. For this purpose, it is important to 
develop schemes that can calculate the fully nonlinear 21cm back- 
ground. 

Given the rapid progress of observations (e.g., GMRT has 
placed an upper bound on the 21cm power spectrum at 2 « 8.6 
in their first result release llPaciga etal.ll20Illl . and MWA and LO- 
FAR are close to their data collection stage), we urgently need a 
thorough understanding of how peculiar velocities enter into pre- 
dictions of the 2Icm signal in observer redshift-space from results 
of modelling or simulations in real space. 



^ http : //www ■ skate lescope . or gi 
^ http://en.wikipedia.org/wiki/Fast_Fourier_Transfomi_Telescope, 
formerly termed Fast Fourier Transform Telescope. 



* I Shaw & LewisI (200 8t) presented a nonlinear analysis of the redshift- 
space distortion. However, they assumed that the 21cm brightness temper- 
ature fluctuati ons are Gaussian, whi ch may not be valid for the EOR (see, 
e.g. Fig. 14 of lMellema et alj2006bh . 
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Along these lines, iMellema et all ( l2006bl . their Figs. 4, 9 and 
10) were the first to consider the effect of peculiar velocities on 
the 21cm brightness temperature fluctuations in observer redshift- 
space when making spectra and maps along the line of sight (LOS). 
They found significant differences between maps and spectra of 
brightness temperature with and without the effects of peculiar ve- 
locities. However, they did not account for this effect when calcu- 
lating statistical properties such as the power spectra of the bright- 
ness temperature fluctuations, nor did they explai n in detail how the 
effects of peculiar velocities were implemented. iLidz et 2i\ ( |2007|) 
claimed to compute the "full redshift-space" 21cm power spec- 
trum, but gave no details of their calculation. Thomas et al. 
claimed to include the effects of peculiar velocities in making 21cm 
maps with their ID radiative transfer simulation, without present- 
ing any details on how these were calculated nor any analysis of the 
effects on statistical quantities. 

The 21cm brightness temperature can diverge in the overdense 
regions of the IGM when corrected for peculiar velocity in the 
optically-thin approximation, because the nonlinear velocity gra- 
dient may cancel the Hubble flow in these regions. In this paper 
we shall investigate the origin of this div ergence and how th is un- 
physical effect can be avoided. Recently, iSantos et al.l ( I2OI0I) pro- 
posed an approximate scheme to circumvent this divergence when 
computing the 21cm power spectru m in semi-numerical models 
of the evolving IGM, also adopted bv lMesinger. Furlanetto & CenI 

l|). In this scheme a numerical cap on the value of the velocity 
gradient is imposed. The accuracy of their approximation, referred 
to henceforth as the "Vu-limited" prescription, has not yet been 
determined. We shall investigate this below. 

Our paper is the first in a series which sets out to build a 
solid and self-consistent computational scheme to predict \he. fully 
nonlinear 21cm background accurately in observer redshift space, 
given density, velocity and ionization fraction information in real 
space. This paper will focus on the methodology for incorporating 
the effects of peculiar velocit y in a nonlinear way. We leave the 
second paper of this series (iShapiro et al.ll20ll]) to focus on the ad- 
ditional nonlinear effects of inhomogeneous reionization coupled 
to peculiar velocity and to test the validity of using the anisotropy 
of the 21cm background fluctuations to separate the astrophysi- 
cal effects of reionization from these of the background cosmol- 
ogy. Some of ou r results were previously summarized by us in 
iMaoetaLlJioiol) . 

This paper is organized as follows. In §|2l we will demonstrate 
how important the effects of peculiar velocity are by comparing the 
angle- averaged power spectra P{k) of brightness temperature fluc- 
tuations when peculiar velocity is neglected, calculated from reion- 
ization simulations, with an approximate scheme that takes pecu- 
liar velocity into account, motivated by linear theory. We then clar- 
ify our terminology in § |3] In § [4] we use a heuristic derivation to 
present a simple picture of redshift-space distortions of the 21cm 
background in the limit of low optical depth and high spin tem- 
perature, and clarify the similarities and differences with galaxy 
redshift surveys. To properly take into account peculiar velocity, 
including the effects of finite spin temperature and optical depth 
and the distinction between thermal- and velocity-broadening of 
the line profile, we present in § |5]the 21cm brightness tempera- 
ture derived from the equation of transfer in an expanding universe. 
We then derive the 21cm power spectrum as measured in redshift- 
space in a hierarchy of approximations, from the exact nonlinear 
power spectrum with finite optical depth to the linear theory in the 
limit of low optical depth. Since the standard linear theory formula 
for 21cm redshift-space distortion (e.g. Barkana & Loeb..2005i) as- 



sumes that all departures from the cosmic mean values (matter den- 
sity, peculiar velocity, and ionization fraction) are of linear ampli- 
tude, while ionization fluctuations are not small for scales compa- 
rable to the size of the H II region, we present here an improved 
version which takes account of ionization fluctuations to higher or- 
der. In § [6] we propose two computational schemes, one based on 
particle data and one based on grid data. We test and compare the 
accuracy and efficiency of these two schemes. In §|7l we investigate 
the accuracy of the optically-thin approximation with regard to the 
21cm power spect rum. In § [H we t est the accuracy of the linear 
theory formula of ISarkana & Loebl (l2005h for redshift-space dis- 
tortion, widely employed to predict the 21cm power spectrum, and 
the new quasi-linear /ik -decomposition presented in §[5] In §|9] we 
discuss the origin of the divergence of the brightness temperature 
found in previous works and how it can be avoided. We also com- 
pare the results of the "Vn-limited" prescription for dealing with 
this problem to the results from our new schemes. We conclude in 
ij llOl We include some technical details of post-processing massive 
numerical particle data in Appendix IaI 



2 HOW IMPORTANT IS PECULIAR VELOCITY? 

Before developing our methodologies we will illustrate the effects 
of peculiar velocities on the 21cm power spectrum as measured in 
redshift-space, to show their importance. For this purpose, we focus 
on the limiting case in which the spin temperature greatly exceeds 
the CMB temperature and the optical depth is small, so we can 
write the differential brightness temperature, 5Ti, = Tt — Tcmb, 
as follows: 



(1) 



where the pre-factor STt is the cosmic mean value in this limit, to 
be defined in equation l l35t . Here Zcos is the cosmological redshift, 
r is the comoving real-space coordinates, and 5p^^ and 5^^ are 
the fluctuations of neutral and total hydrogen density in real space, 
respectively; i.e. pm = Ph xm, and^p^^j = Sp^+Sx^i+^PH ^^hi' 
where xm is the neutral hydrogen fraction. Also, we define the 
quantity 



1 + Zcoa dvu 



dr\ 



(r), 



(2) 



the gradient of the proper radial peculiar velocity along the LOS, 
normalized by the conformal Hubble constant H/{1 + jZcos). 
The power spectrum of brightness temperature fluctuations in ob- 
server redshift-space can then be written as ^5T^* (k)5ri,(k')^ = 

(27r)^Pl?,(k)5'^)(k - k'), where 5ft(k) is the Fourier transform 
of 5Tb. Hereafter Px.x is the auto-power spectrum of the field x, 
and Px,y is the cross-power spectrum of the fields x and y. 

In Figure [T] we present slices for three versions of the three- 
dimensional power spectrum, the first being the one without includ- 
ing any effects of peculiar velocities (hereafter dubbed the "uncor- 
rected for peculiar velocity", or UPV, scheme), given by 



3UPV,3D 



(k) = 5n(zc, 



k). 



(3) 



The second version is calculated according to the "quasi-linear 
/Xk-decomp osition scheme" (a generalization of linear theory in 
iBarkana & Loebi.2005i : but see the exact definition and derivation 
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kx (h/Mpc) 



Figure 1. 3D power spectra A^(k) = A;^P2i (k)/27r^ (in units of mK^) of 21cm brightness temperature fluctuations. Tlie panels show a slice through the 
kx-ky plane, with the LOS along the x-axis, calculated from our numerical simulation at the 50% ionized epoch (z = 9.457). Top left: UPV scheme; top 
right: quasi-linear ^^"decomposition scheme; bottom: the fully nonlinear PPM-RRM scheme. 



in §[53]beIow), 

On large scales, according to linear theory, the second and fourth 
moments of the /ik-decomposition in equation Q come from the 
cross-coiTelation of the peculiar velocity gradient with neutral hy- 
drogen density fluctuations and the auto-correlation of the peculiar 
velocity gradient, respectively. Here /Xk = fe||/|k| where fcy is the 
LOS component of k. The moments in the RHS of equation ^ 
are angle-averaged in a spherical k-space shell with k = |k|, i.e., 



Ps^ (k) = (Ps^ S'- fk)), etc. Note that in equation (HI, 

PHI' Pm PHI' PHI ^ " 1 " 

the quasi-linear /ik-decomposition power spectrum can be com- 
puted directly from the real-space data, avoiding the need to specify 
a computational scheme for calculating the redshift-space-distorted 
21cm signal data cube. 

The third version of the three-dimensional power spectrum 
of 21cm brightness temperature fluctuations shown in Figure [T]is 
calculated using a numerical scheme that finds the fully nonlin- 
ear redshift-space-distorted 21cm brightness temperature signal as 
a function of position and frequency (the "PPM-RRM" scheme, see 
§[6j2jTJ. This last version of PatO^) will be derived in the sections 
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Figure 2. Ratio of 21cm redshift-space-distoited power spectram in tlie 
quasi-linear /^^-decomposition scheme and 21cm power spectram in the 
UPV scheme as a function of comoving wavenumber k, at a series of red- 
shift z and mass-averaged ionization fraction Xi jvf ■ Arrows indicate the 
direction of the evolution of the curves at low k as reionization proceeds. 
Starting from the curve at Xi^M = 0.002 (dark blue, short dash - long 
dash) near the ratio = 1.87 limit, the ratio at low k moves up through 
the curves, in sequence, at M = 0.009 (purple, short dash - long dash), 
Xi,M ~ 0.05 (orange, short dash - long dash), m 0.1 (dark red, short 
dash - long dash), flips the direction at x^ j^i « 0.2 (magenta, dot - long 
dash), then moves down through the curves at Xi jvf ~ 0.3 (cyan, dot - short 
dash), flips the direction again at x^ m ^ 0.4 (blue, long dash), moves up 
through the curves at Xi^M ~ 0.5 (green, short dash), Xi^M ~ 0.75 (red, 
dot), and approaches the curve at x^ m ^ 0.9 (black, solid) near the ratio 
= 1 limit. 



which follow, based on the results of numerical reionization simu- 
lations. 

The simulation data used for Figure [T| are taken from a ra- 
diative transfer (RT) simulation of a 114 h^^ Mpc box with 256^ 



RT resolution (more fully presented in ij |6.3b . For the UPV scheme 
(top left), the power spectrum is seen to be numerically fluctuating 
in equal-|k| shells, but otherwise to be isotropic in the sense that it 
does not show any directional preference. For the quasi-linear /.ik- 
decomposition scheme (top right), the power spectrum is perfectly 
distorted along the LOS direction, i.e. elongated for small k and 
squeezed for large k. For the PPM-RRM scheme (bottom center), 
it is hard to see the distortion for the small- A: modes due to the small 
number of modes, but the compressed nature of the large-fc modes 
is clearly visible, albeit with some numerical noise. Clearly, pecu- 
liar velocities introduce noticeable anisotropies in the 21cm power 
spectra. 

To make a more quantitative comparison, we compute the 



.qlin,3D 
AT 



angle-averaged power spectrum P^y''°'^°(fc) 
for the quasi-linear /ik -decomposition scheme, 

3 Pn PHI ^ ' 5 PH ' PH ^ ' 



(k)) 



(5) 



as a function of k = |k|, and the same for the UPV scheme, 
ST^{zcoa)Ps^ (fe). Figure |2] shows the 

" ^ ' Out ' Out ^ ' ^ 



ratio of these two, P^'^""'^°(fe)/P^y^'^°(fc), for ten different 
phases of reionization. Two limiting cases are obvious: for the 
early phases of reionization, the ratio approaches an almost con- 
stant value of 1.87; for the late phases the ratio tends to 1.0. These 
limits hold best at low k. They can be understood as follows. At 
early times, the neutral fraction fluctuations 5^.^-^ are negligible, 
i.e. the neutral hydrogen density traces the total hydrogen density 
almost exactly, and, therefore, the 21cm power spectrum in the 
quasi-linear /ik-decomposition scheme differs from P^^^'^° ~ 

-2 



by a factor of 1 + ^ + f = 1-87. At late times. 



neutral fraction fluctuations dominate over density fluctuations Q , 
so its auto-power Pgr gr becomes the dominant term in both 
versions of the power spectra, making their ratio approach unity. 
As pointed out above, the density fluctuation terms in equation 
l|4) reflect the redshift-space distortion caused by peculiar veloc- 
ity. Hence, the effect of peculiar velocity on the power spectrum 
of 21cm brightness temperature fluctuations becomes subd ominant 
towards th e end of reionization, as note d also by McQuin n et al.| 
1 I2OO6I) and lMesinger & Furlanettol ( |2007|) . 

Between these two limits. Figure |2] shows that as reionization 
proceeds the ratio evolves rather nonlinearly, changing both ampli- 
tude and shape non-monotonically. Reionization proceeds "inside- 
out" in our simulation, i.e. overdense regions ionize earlier than 
underdense regions, so the cross-power Pgr s'- between den- 

PH' ^HI 

sity fluctuation and neutral fraction fluctuation is negative at large 
scales. Shortly after the onset of reionization (xi < 0.2), the to- 
tal density power spectrum Pgr gr still dominates over the other 

terms, but the cross-power Pgr 51- also contributes significantly 

ph' ^hi 

and is the next most important term, so P^'^''"'^'^/P^^^'^'^ 



(1.87P5. ,6. + 2.67Psr ,sr )/iPsr 

^ PH' PH PH ' ^HI " ^ PH' PH 



AT 
+ 2P«- ,Sr 



1.87 - 1.07(P5r- gr 

the neutral fraction 



/Pgr gr ), moving the ratio up since 

' PH ' PH ^ ^ 



uctuations increase as reionization pro- 
ceeds. When reionization reaches the midway point (xi > 
0.4) and large ionized bubbles have formed, the neutral frac- 
tion auto-power Pgr gr starts dominating over other pow- 
ers and the cross-power Pgr gr becomes subleadins, so 

PH' ^HI ° 

pa,qlin,lD ;pUPV,lD 
AT /-'AT 

m)iPs^.^,Sr_/Pgr 



l+(2/3)(P,. 



1 + 



PH '"^HI '^"^'^PHI ''^PHI ^ 

gr ). Since the cross power Pgr gr 

PH'"^HI' "==HI' ^Hl' PH ' ^HI 

is negative, the ratio is less than unity. As reionization proceeds 
towards its final stages, the neutral fraction fluctuations continue 
to grow, pushing the ratio closer and closer to unity. Between 
Xi ~ 0.2 and 0.4, the competition between neutral fraction fluc- 
tuations and density fluctuations makes the ratio at large scales first 
turn around at a large value ~ 4 — 5 (xi « 0.2), then move all the 
way down to less than 1 (Xi ~ 0.4), then turn around again and 
begin to approach the limit of 1 . 

These comparisons illustrate the nontrivial effects when 
applying redshift space distortions using the quasi-linear /ik- 
decomposition scheme. However, the fully nonlinearly distorted 
21cm power spectrum may well show a more complicated behav- 
ior, which has not been previously explored. In order to calculate 



^ The variance in neutral fraction can be estimated as (Axm)^ = 
({^^HI ~ ^ijji)^) ~ 5;hi(1 — Shi), so the rms neutral fraction ^wrtwa- 



tion is 5™= = AxHi/a;Hi 



xjii)/xjii. Thus the neutral frac- 
tion fluctuations grow as reionization proceeds, even though the variance in 
neutral fraction decreases near the end of reionization. 
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the fully-nonlinear redshift-space-distorted 21cm power spectrum, 
we need a robust scheme to compute it from simulation results. The 
aim of this paper is to develop such a scheme, taking into account 
all peculiar velocity effects. In a subsequent paper we will use this 
scheme to study the nonlinear distortion in the 21cm power spec- 
trum and test the validity of quasi-linear ^k-decomposition scheme 
upon which the 21cm cosmology is based. 



3 TERMINOLOGY 

Before we proceed to the main content, we summarize in this sec- 
tion our terminology which otherwise may be confusing. 

3.1 Reference Frames 

We distinguish between different reference frames. These are 

• Emitter space: the local rest-frame of the emitting atoms. 

• FRW space: the cosmic reference frame in which space is 
uniformly expanding, as described by the Friedmann-Robertson- 
Walker metric. 

The emitter space and the FRW space are related by the lo- 
cal Lorentz transformation at the position of emitting atoms, and 
the relative motion of these two frames is the peculiar velocity of 
atoms. 

From the observer's point of view, the coordinates of source 
{t, r) in FRW space can be relabeled by tarrivai (arrival time of ra- 
diation emitted at time t by source located at comoving location r), 
Zcos (cosmological redshift experienced by photons from time t of 
their emission to the time tarrivai at which they reach the observer), 
and (angular coordinates on the sky). 

• Observer real space: for fixed tarrivai = tprcscnt (present 
time), the observer can reconstruct a part of the FRW space theo- 
retically ~ those (t, r) on the light-cone that can be determined by 
Zcos. In particular, r = ''(zcos)jip„,„„t = Jq""' cdz'/H(z'). We 
call this the observer real space. In the rest of this paper, quantities 
measured in real space are superscripted with r, so for example 

is a number density in real space. 

• Observer redshift space: in practice, observers can only mea- 
sure the observed redshift of radiation, since the wavelength is red- 
shifted both cosmologically and by the Doppler shift associated 
with peculiar velocity, z/obs = 2^0/(1 + ^obs) and 1 -I- Zobs = 
( 1 -I- Zcos )(1 — ^)~^- Observers can set up a "distorted" comoving 
coordinate system, known as observer redshift space, in which the 
position of the emitter is the apparent comoving position if the red- 
shift is interpreted as only cosmological, i.e. s = 7"(zobs) Itprcacnt = 
f^°^' cdz' /H{z'), which shifts the real comoving coordinate r 
along the LOS (f) to 

I (-^ ^ Zobs) /i \ ' //-x 

s^r+—— —vii(t,r)r. (6) 

-n (Zobsj 

Note that the transformation between observer real and redshift 
spaces is not covariant (in a general relativistic sense) or even 
Galilean invariant, since it does not preserve spatial intervals at 
fixed time. In the rest of this paper, quantities measured in observer 
redshift space are superscripted with s, so for example, n" is a num- 
ber density in redshift space. 

3.2 3D Mapping Distortion 

One can distinguish between several types of distortions, namely 



• Apparent location distortion in redshift-space: when the 
observed frequency of a spectral line from a distant source is used 
to locate the source along the LOS, the answer depends upon solv- 
ing equation which requires a knowledge of the LOS peculiar 
velocity of the source at the time of emission. The term "redshift- 
space distortion" usually refers to the error one makes in locating 
the source by assuming the peculiar velocity to be zero. 

• Brightness temperature distortion in real-space: Radiative 
transfer effects can result in a modification of the observed 21cm 
brightness temperature due to gradients in the velocity field, as 
shown in § 15.11 This effect is independent of the adoption of ei- 
ther real- or redshift-space. In other words, even if an observer 
could construct a 3D mapping of 21cm brightness temperature in 
observer real space by knowing the peculiar velocities along the 
LOS, gradients in the peculiar velocity field can still modify the 
magnitude of brightness temperature. 

• 21cm redshift-space distortion: This is the combination of 
the previous two distortions, namely the apparent location distor- 
tion in redshift-space and the brightness temperature distortion in 
real-space. The observed 21cm signal is modified by the presence 
of peculiar velocities according to this combination. 



3.3 Power Spectra 

Power spectra can be calculated in different dimensions in k-space 
and with different methods for applying the effects of peculiar ve- 
locities. We use the following terminology: 

• 3D power spectrum P3D (k) : The power spectrum in three- 
dimensional k-space. 

• ID power spectrum PiD(fc): The power spectrum in one- 
dimensional |k|-space (or simply fc-space), obtained by averaging 
the 3D power spectrum over modes in spherical shells in k-space: 

PlD(fc) = (P3D(k))^,^^ll Withfc=|kl. 

• 21cm power spectrum: An abbreviation of "power spectrum 
of 21cm brightness temperature fluctuations". 

• 21cm redshift-space-distorted power spectrum: The 21cm 
power spectrum in observer redshift space, i.e. taking the 21cm 
redshift-space distortion into account. 

• 21cm real-space power spectrum: The 21cm power spec- 
trum evaluated with velocity gradient coiTections and yet in real 
space, i.e. the power spectrum which results from the Fourier trans- 
form of the scalar field corresponding to the true (i.e. peculiar- 
velocity-corrected) 21cm brightness temperature at each point in 
real-space at a single cosmic time. This power spectrum so-defined 
is not the power spectrum of the observed 21cm brightness temper- 
ature field evaluated in redshift-space in which each plane trans- 
verse to the line-of-sight corresponds to a single observed fre- 
quency. Instead, this "real-space power spectrum" represents the 
brightness temperature at different observed frequencies for differ- 
ent locations in real space, as a result of Doppler shifts caused by 
peculiar velocity. 

• 21cm UPV power spectrum: The 21cm power spectrum 
evaluated without any velocity gradient corrections and in real 
space, i.e. taking into account neither the brightness temperature 
distortion in real-space nor the apparent location distortion in 
redshift-space; "UPV" stands for "uncorrected for peculiar veloc- 
ity". 

• 21cm quasi-linear -decomposition power spectrum: an 

abbreviation of 21cm power spectrum calculated with the "quasi- 
linear /ik-decomposition scheme" (see ij |3.4t . 
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3.4 Computational Schemes 

Below we develop different schemes for applying the effects of pe- 
culiar velocities. We summarize these here. 

• Linear theory: A scheme to compute 21cm power spectrum 
in redshift-space, where all fields, de nsity, velocity and ioniz ation 
fraction, are linearized; introduced bv lBarkana & Loel? (2005). 

• Quasi-linear /.1,^-decomposition scheme: A scheme to com- 
pute 21cm power spectrum in redshift-space, assuming the density 
and velocity fields to be linear, but without constraints on the ion- 
ization fraction field; introduced in § 15.3.31 

• PPM-RRM scheme ("Particle-to-Particle-to-Mesh Real-to- 
Redshift-Space-Mapping"): A particle-based numerical scheme to 
construct the 21cm data cube in observer redshift space, using the 
real- to redshift-space re-mapping of density, velocity and ioniza- 
tion fraction data. Introduced in § 16.2.11 

• MM-RRM scheme ("Mesh-to-Mesh Real-to-Redshift- 
Space-Mapping"): Same as the PPM-RRM scheme, but grid- 
based; introduced in § 16.2.21 

• DEMRF scheme ("Direct Evaluation by Multiple Real-space 
FFTs"): A scheme to compute the 21cm power spectrum in 
redshift-space by a direct integration technique; introduced in 
§1133] 



4 21CM REDSHIFT SPACE DISTORTION: OPTICALLY 
THIN AND HIGH Ts LIMIT 

In this section we consider the simplest scenario, namely in the 
limit of small optical depth and high spin temperature Ts ^ Tcmb , 
and show that in this limit, peculiar velocities affect the 21cm 
brightness temperature in an analogous way to the redshift-space 
distortion in galaxy redshift surveys. 

Recall that galaxy redshift surveys can distinguish individual 
galaxies. In other words, galaxies can be counted directly. Peculiar 
velocities move galaxies to their apparent locations, thereby affect- 
ing the number density of galaxies in redshift-space. For 21cm sur- 
veys, however, individual 21cm-line emitters — each neutral hy- 
drogen atom — cannot be resolved and the H I number density 
can only be inferred from the observed brightness temperature of 
21cm emission. This fundamental difference from galaxy redshift 
surveys implies that radiative transfer effects associated with pe- 
culiar velocities must be taken into account when calculating the 
redshift-space distortion of the 21cm background. 

In the optically thin limit, the emission from each individ- 
ual H I atom can be regarded as independently transferred along 
the LOS. In the high spin temperature (Ts 3> Tcmb) limit, the 
stimulated emission/absorption is negligible compared to the spon- 
taneous emission. Therefore when both limits apply, each H I 
atom can be thought of as an independently shining 21 cm-line 
source with the intrinsic luminosity = huoAu), where — 
21 cm/c = 1420.4057 MHz, and ^lo = 2.85 x lO^-^^s"^ is the 
Einstein spontaneous emission coefficient of the 21cm transition. 
Then the emissivity at frequency i/pp in the emitter space is 



47r 



(7) 



where n[ ~ (3/4)nHi is the number density of H I atoms in the 
upper hyperfine state in real-space. The function 4'{'^rf) is the line 
profile and satisfies the normalization condition 4'(.^) = 1. 
The radiative transfer equation in FRW space then becomes 



(8) 



where is the comoving specific intensity of a light ray. The ray 
path can be labeled by the proper distance along it, = cdt, 
where t i s the physical t ime. Since the emissivity transforms as 
(see, e.g.. lMihalasll 19781) . > = (1 ~ j^'^ in FRW space. The 
observed specific intensity then is 



47r 



■m<^(i/R,F)d^. 



(9) 



In the idealized case of no thermal broadening, the line profile 
is a (5-function peaked at the transition frequency seen from the 
emitter space, i.e.. 



<^(i^Rf) = 5(l^RF - fo) . 



(10) 



Therefore, the integration picks up the integrand evaluated at the 
location of emission. Using an identity, whose derivation will be 
described in §[5] 



di 



= -ua H(a) 
c 



1 + 



dv\\ 



aH{a) cir|[ 



(11) 



where ry is the comoving LOS distance, we find that the observed 
specific intensity is 



(12) 



1 + 



aH{a) dr| 



Now we consider the distortion of apparent location. The num- 
ber density in redshift-space nl satisfies 



ni 



1 + 



(13) 



aH(a) dru 



since the number of H I atoms are preserved between real- and 
redshift-space, and the volume element in redshift-space is dis- 



torted as SV = SV 1 + 



dvu 



aH(a) drii 



. Therefore we find that 



47r !/gi?(a) 



ni 



(14) 



(We will express eq. [14] in terms of the familiar brightness tem- 
perature in §[5]) This is to say, the simple proportionality relation 
between the specific intensity (or, brightness temperature) and the 
neutral hydrogen density is preserved with and without peculiar ve- 
locities. 

There is a simple explanation for equation ( I14t . In the limit of 
optically thin and high spin temperature, 21cm radiation from each 
neutral atom is emitted and then transferred independently. There- 
fore, the radiative transfer effects of peculiar velocity on a pocket 
of gas is simply equivalent to the simple picture of having all emit- 
ters shine from their apparent locations. Note that this net effect 
combines the peculiar velocity effects on the radiative transfer and 
on the distortion of apparent locations of sources. 

Equation l ll4t establishes that, in this limit, there is an analogy 
between 21cm brightness temperature measurements and galaxy 
number density measurements, in that the neutral hydrogen atom 
number in 21cm tomography corresponds to the galaxy number in 
galaxy surveys. Therefore the 21cm power spectrum should be af- 
fected by peculiar velocities in a form similar to the linear redshift 
space distortion on large scales (first shown by Barkana & Loebl 
l2005h . similar to the galaxy matter power spectrum dKaiserll 19870 . 
In both cases, the effects of peculiar velocities can be thought of 
as the distortion due to displacing sources to their apparent LOS 
locations. 

In the more general case in which optical depth is not small 
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and/or Ts < Tcmb, however, the analogy between the redshift- 
space distortion of the 21cm background signal and that in galaxy 
redshift surveys breaks down. Since galaxy redshift surveys can 
resolve and count discrete galaxies, they do not depend upon mea- 
suring the unresolved intensity of galactic emission to deduce the 
number density of galaxies. For the 21cm background, however, 
we cannot resolve individual sources, (i.e. individual atoms), so we 
must use the specific intensity (or brightness temperature) to infer 
the source density, e.g. in the optically-thin/high Ts limit, accord- 
ing to equation l ll4t above. If the conditions of low optical depth 
and high spin temperature are not satisfied, however, equation l ll4t 
no longer applies. In that case, the luminosity emitted per atom then 
depends upon the unknown spin temperature, and the received in- 
tensity is also no longer linear in the optical depth. In order to inter- 
pret redshift-space-distorted 21cm maps, in general, therefore, we 
cannot simply borrow the analogy of the galaxy redshift surveys. 
We discuss the details of this in §(5] 



5.1.1 The Formal Solution 



5 EFFECTS OF PECULIAR VELOCITY ON THE 
OBSERVED 21CM BACKGROUND 

Given density, velocity, and ionization information in real space, 
peculiar velocities can affect the observed 21cm signal through two 
effects: (1) the observed 21cm brightness temperature can be mod- 
ified by the gradient of radial peculiar velocity of the gas along the 
LOS, and (2) the apparent location of the gas can be shifted from its 
reaZ-space location because of the Doppler shift due to its peculiar 
velocity. We will address the first effect in i) l5.1l and then combine 
both effects to form a self-consistent picture of 21cm redshift space 
distortion in ij|5.2|and ^\5.3\ 



5.1 



The Transfer of 21cm Radiation Tlirough the 
Intergalactic Medium 



The effect of peculiar velocity gradients on obs erved 21c m bright- 
ness t emperature was first addre ssed in iBharad wai. Nat h_& Sethil 

l iooil) 



and 



Bharadwai & All 



subsequently _^ 

Barkana & Loeb (2005). IBharadw ai. Nath & Sethi] ( 1200 ll 



Bharadwai & Alii ( 120041) only explored the simpler limit of high 
spin temperature (Ts ^ Icmb) and optically thin radiative trans- 
fer, and implicitly assumed that the velocity gradient is small so that 

the factor 1/^1-1- ^g^^,) ^Tjj") '-^ linearized. iBarkana & Loebl 
( l2005h attributed this velocity gradient correction to the effect of 
the fixed thermal width of the 21cm scattering cross section, with- 
out showing the details of the derivation. Since we aim to under- 
stand peculiar velocity thoroughly, it is worthwhile to re-derive this 
effect from first principle, i.e., solving the radiative transfer equa- 
tion, and keeping all contributions of peculiar velocity to linear or- 
der v/c. In this section we show that it is the peculiar velocity of 
the hulk motion, not the thermal broadening, that is responsible for 
making its correction in 21cm brightness temperature. We check 
the validity of the optically thin approximation and show that it 
can break down in certain conditions, although it is mostly valid in 
the IGM. We find also that, in addition to the well-known velocity 
gradient correction, the contribution of spin temperature to 21cm 
brightness temperature can be modified by a term of order 0{v/c). 



Consider a light ray with comoving specific intensity Iv passing 
through a gas element. In an expanding universe, in which v oc 
l/g, the radiative tr ansfer e quation reads (Gnedin & O striker. 19971 : 
Iwise & AbefcoilUzhangTHui & Haimanll2007i) 



div h 
cadrj a 



VJ,. - 



H{a) dh 
c dlni/ 



(15) 



where Ii, is a function of conformal time rj, comoving coordinates 
r, frequency u and direction h. Here a is the cosmic scale factor 
and H{a) is the Hubble constant at a. The ray path can be labeled 
by the proper distance along it, — c dt, where t is the physical 
time. The radiative transfer equation can be rewritten in terms of 
the Lagrangian total derivative 



dh 
d^ 



I^l/Iu + jl/ ■ 



(16) 



Here and are the absorption coefficient and the comoving 
spontaneous emission coefficient at the frequency in FRW space, 
respectively. 

We label i^obs to be the frequency observed today, v' — 
i^ohs/a. the frequency at some proper distance 5' along the ray path 
in FRW space, and u^p — ^)^^ the frequency in emit- 
ter space, where «|| is the radial proper peculiar velocity of the gas. 
Hereafter, the subscript or superscript "RF" stands for "rest-frame". 

By defining the optical depth forward along the ray path as 



dr^, = K^/ d^' , 



(17) 



the radiative transfer equation has the formal solution for the spe- 
cific intensity observed today at frequency i^obs 



r _ rCMB - 



(18) 

Here we assume that the ray has the same comoving specific inten- 
sity as the CMB (Z^?"^™) when the ray was on the far side of the 
gas element from the observer. S^i{^') — j„' / Hi,' is the comov- 
ing source function at the frequency v' seen in FRW space at the 
proper distance ^' on the ray path. Ti,^^^^ is the integrated optical 
depth through the gas. 



5.1.2 Optical Depth 

In emitter space, the absorption coefficient is 

= -huo{noBoi — ni_Bio)0(t'Rp) , 
c 

where Bqi and Bio are the Einstein probability coefficients for 
induced upward and downward transitions, respectively, between 
the lower state with density no and higher state with density ni. 
The spin temperature is defined to be the excitation temperature 
between the hyperfine states, i.e. 
m 
no 



(19) 



Si -T^/Ts 

— e 



3e 



go 



(20) 



where go — 1 and gi — 3 are the statistical weights. = 
hi^o/ks = 0.068 K is the temperature corresponding to the rest- 
frame frequency vq. For 21cm transitions, all astrophysical appli- 
cations satisfy Ts 3> T*, so no = nHi/4 where nni is the proper 



It is sometimes customary to use the proper specific intensity III 



(p) 



which is related to the comoving specific intensity by = /, 



(P)„3 
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number density of neutral hydrogen. 01t is straightforward to show 
that 



RF _ 3c^ylioT^nHi<^(^'R.F) 



(21) 



using the identities giBio — goBoi = c^giAio/SnhuQ. 

Now we transform our calculation to FRW space. The absorp- 
tion coefficient transforms as (see, e.g.. .Mihalas, , 1 97 . so in 
FRW space 



RF/ / / /-. RF '^ll N- 



(22) 



Therefore the optical depth is 



'^(i'rf) d£_' . (23) 



For the 21cm line transition, in the idealized case of no ther- 
mal broadening, the line profile is (/'(t'Rp) = 5{^'rf ^ ^o)- Inte- 
grating a (^-function picks up the integrand evaluated at the peak 
which physically corresponds to the location on the ray path where 
the transition actually takes place, its proper distance labeled as 
(hereafter in this section, the subscript r stands for "radiation"). We 
assume that each ray with a given observed frequency only experi- 
ences one 21cm transition event along the ray path. (We discuss the 
multi-transition case in i; 15. 1.71 ) The line profile can be rewritten as 



\dv'^^/d£_'\ 



(24) 



Here we assume the non-singular case, i.e., (di^Rp/i^^Oj 7^ 0. 
(We discuss the singular case in § 15.1.51 ) We use the relation 

fRF = i'obsa~^(l r)"^ to take the derivative di/^p/d^', and 

then evaluate it at where j/pp = z^o- It is straightforward to show 
that 



flr c drii 



(25) 



where V\\ — ar\\H{a) + v\\ is the proper velocity along the LOS 
and 



dm 



aH{a) + 



dr« 



(26) 



with ry the comoving LOS distance. Therefore the optical depth is 

_ 3c^AloT-»arWHl(^r) 



327ri.gr,(Cr) |9Vj|/ar||l (1 



5.1.3 Observed Brightness Temperature 



(27) 



Now we simplify the formal solution of radiative trans- 
fer equation. Since dr^/ oc — ^r)d^', the inte- 
gral in equation dlSt takes non-zero contribution only from 
^' = Cr, therefore 5',, (C')e"*""°'- = 
S.'{&) Jo"""' e-(-'^ob=-<')dr:, = S,,{^r){l - e'-'^obO.H 

tt Note that strictly speaking, nni is the number density in emitter space. 
But the number densities in emitter space and in FRW space only differ 
in the relativistic limit, i.e., n^^/n'^™ = dVcos/dVjip = dtjip/dt = 
1 / ^1 — t)^/c^, so we can ignore the difference to linear order v/c. 



The factor e 



is a step function at 



more rigorously, the integral yields Jq"°^' S^:{£,')e '^"obo '^v'^dr^, = 
'S'l/' (?r)T'i^obB [-t ~ (1 ~ e^^^oba )r)(0)]. The unit step function ri(x) at 



In emitter space, S^f = 2kBVQTs{£,r) / , i.e. the Planck 
function evaluated with the spin temperature Ts at The source 
function transforms as (see, e.g.. Mihalas,, 1978i) . so the comov- 
ing source function in FRW space is 



(28) 



where accounts for the comoving factor. 

Suppose the ray has the frequency Vp — Uohs/op with some 
scale factor < ar, (i.e., when it is on the far side of the gas 
element from the observer,) and is in equilibrium with the CMB of 
temperature rcMB,p = Tcmb,o /cip. In the absence of intervening 
atoms, the comoving specific intensity observed today would be 

I^^^ ~ a^2kBVpTcMB,p/c^ = 2fc_gt'obs2cMB,o/c^. 

The 21cm brightness temperature at the observed frequency 
I'obs is defined by 



ThiVoh 



From equation ( I18t it is straightforward to show that 



-)(1 



(29) 



(30) 



The 21cm line is generally optically thin to the IGM, i.e. t^^^^^ <^ 
1. (We discuss the validity of this approximation in 15. 1.61 ) In this 
limit, the differential brightness temperature is 



STb{vohs 



Tb{l'ohs) — rcMB,o 

3c^AionnHi(r)ar 
32^r.3//(aO|l + (aJ/)-i^(r)| 

rCMB(tSr) 



(31) 

i{ar)] , (32) 



(33) 



where r is the real-space location of 21cm transition corresponding 
to the proper distance on the ray path. Tcmb (fflr) = Tcmb.o / cir 
is the CMB temperature at the time of 21cm transition. Here we 
define the ejfective spin temperature 



rf(r) 



r.(r) 



1 



^11 (r) 



(34) 



Note that equation ([34} only infers that the spin temperature man- 
ifests itself to 21cm brightness temperature and optical depth in 
a manner modified by the peculiar velocity, but this effect does 
not modify the level population of hydrogen hyperfine states, 
nor the spin temperature. The level population can in fact be 
modified by pec uliar velocity through an effect pointed out by 
IChuzhov & Shap iro (2006). This is however a different effect from 
the one in equation ([34} which is based on a given spin temperature. 

Equation (I33t i s in ag reement with the well-known equation 
in iBarkana & Loebl ( |2005|) except for the appearance of effective 
spin temperature T^^ . However, this modification is actually not 
important for two reasons. First, it is of order 0(v/c) and the bulk 
motion of gas is mostly non-relativistic. Second, many research pa- 
pers focus on the epoch during reionization when Tcmb/Ts <Si 1, 



i: = is undefined in general, but using an identity intrinsic in this prob- 
lem 1 - e-^-oba = /^-o^- g-t-'ob. = r.^^Jl _ (1 - 
e~'^'^ob!i )ri{p)\, we can regulate r;(0) and obtain the same result. 
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when the spin temperature has a neghgible effect on the brightness 
temperature. 

For convenience, we define the me aiHbri ghtness temperature 
in the hmit Ts 3> Tcmb as 



23.i 



3c^y4ior*nHi(zcos) 



^COS 



0.02 



0.15 l + 2c 

Om/i^ 10 



(2:cos)mK, (35) 



where the cosmological redshift is defined by 1 + Zcos = 1/a,- 
nmizcos) is the mean neutral hydrogen number density at z. 
xm,-m{zcos) is the mean mass-weighted neutral fraction at z. 
Then we can rewrite equation ( 133) in terms of fluctuations. 



cos, 
COS . 



l + 5n 



|l + 5s,„(r)| 



1 - 



7cMB(ar) 



TS«{v) 



(36) 



where (5pjjj (r) = [7iHi(r) — nHi(zcos)]/nHi(2cos) is the neutral hy- 
drogen density fluctuation, and 5a^i,(r) is defined in equation 



5.1.4 Line Profile Revisited: Velocity vs. Thermal Broadening 

In general, the line profile can include thermal broadening, as well 
as velocity broadening due to bulk motion. The velocity broadening 
is naturally included by taking the 5-function-shaped line profile 
peaked at the rest-frame frequency which is shifted both cosmo- 
logically and by Doppler effect. Our calculation (eqs. |27|and|36|l 
shows that the velocity gradient correction is due t o the bulk motion 
of neu tral atoms. However, in their original paper. Isarkana & Loebl 
( l2005h explained the inclusion of the velocity gradient compactly, 
without showing details, as "The velocity gradient term arises be- 
cause the 21 cm scattering cross section has a fixed thermal width, 
which translates through the redshift factor (1 + Hr/c) to a fixed in- 
terval in velocity". This seems to mean that the thermal broadening 
is responsible for the velocity gradient correction. In this subsec- 
tion, we will clarify that in the non-singular case, the contribution 
of thermal broadening is always subdominant to the velocity broad- 
ening. 

Basically, the thermal velocity of hydrogen atoms can con- 
tribute an additional Doppler shift of the line frequency. For a given 
j^obs, neutral atoms can in principle have a possibility, given by the 
Maxwellian distribution, of seeing the radiation in the 21cm rest- 
frame frequency vo, even if i/^p (with Doppler shifted due to bulk 
motion) 7^ uo . This is described by the Gaussian line profile, re- 
placing equation l llOt . 



0(;^rf) = 



exp 



where 



c 



2kBTk 



(37) 



(38) 



is the thermal Doppler shift corresponding to a gas kinetic temper- 
ature Tfc. 

In the non-singular case, i.e. {dv^-p /d^')^ 7^ 0, we can 
change the integration variable in equation ( I23t . 



(39) 



This is not a volume-weighted mean, but essentially a mean in the 
redshift-space. See footnote l24l 



and rewrite the optical depth with thermal broadening as 



(40) 



where T(5') is the function in equation l l27t with replaced by ^' 
corresponding to J^rf, so by definition T(^r) = ^i^^^ is the optical 
depth without thermal broadening for the observed frequency i^obs. 

Suppose the rest-frame frequency finds i^^j^F ~ '^0 at ^' = ^r. 
The thermal width is small compared to uq, since Afth/j^o ~ 10~^ 
if Tfc ~ 10^ K. Therefore the integrand is nonzero only near u^p — 
1^0. So we can Taylor expand the integrand at vo to sub-leading 
since the first order oc di^^F('^RF ~ 

= 0. It is straightforward to show that the 



order in 0{v'p^p — vq)^, 
-)exp[-%i^] 



result is 



T NT r, , A T 1 

where the fractional correction due to thermal broadening is 



(41) 



d^l 



dvL 



A 2 



Afth 



. 10-^ 
(42) 

Here we assume the gas temperature is about 10* K, which is close 
to the maximum temperature attainable by neutral hydrogen before 
collisional ionization becomes important. Therefore, the contribu- 
tion of thermal broadening is always negligible compared to the 
bulk motion. 



5.1.5 Obserx'ed Brightness Temperature: Optically Thick Limit 

Our results for optical depth and brightness temperature seem to 
diverge for 5g^,^ = —1 (see eqs. |27|and|36t. We discuss this sin- 
gularity behavior in this subsection, and find that the divergence 
in optical depth can be relaxed by including thermal broadening, 
and the divergence in brightness temperature can be removed by 
dropping the optically thin approximation. 

We should first note that the singularity at 5a^„ = — 1 cor- 
responds to {dv^p/dS,')^ — 0. In this case, the regular changing 
variable technique (eq.|39t is invalid. Instead, one should Taylor ex- 
pand J'rf(C') iisar to second order, i/^p — uo + ^ ^r)^' 



where /? 



, and find that 



de' = sgn(C'-6-)sgn(/3) 



Then the optical depth becomes 

a(/3)-cx) 

di^RF sgn(/?) 



= 2 



Sc^AioT.nHi 



(43) 



(44) 



If there is no thermal broadening, the line profile is a 5- 
function peaked at v'p^p — vo, and the optical depth is still divergent 
due to the 1 / \/ ^^rf — 1^0 factor However, thermal broadening can 
remove this divergence. To see this, we can move the ^'-dependent 
factors (nm, Ts, and ny) out of the integral, evaluated at ^r, since 
the evaluation is concentrated near ^r. When applying the thermal 
line profile (eq. |37}, we find that the optical depth in the singular 
case {Sg^v = — 1) becomes 

Sc^AioT^nmier) 1.446 



ii(«.) 



) Vm 



(45) 



t^th 
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Here the factor 1.446 is an approximation of \/2/it x 2r(5/4). 
Since t^^^^^ oc l/v^Ai^, the optical depth can be large when Sg^^ 
is close to —1. We evaluate /3 when 5g^^ = — 1, 



/3 = - 



uoH{ar 



2 + 



arH'{ar) c 
H{ar) {arH{ar)Y dr?. 



(46) 



where H'{a) = dH/da. 

There are two astrophysical cases that can generate Sg^v ~ 
— 1. One is the virialized halo and the other is the spherical col- 
lapse at the turn-around point (pre-virialization). In both cases, the 
proper velocity is V|| — ar<^\H[a) + =0 (seen from the halo 
center), so dv\\/dr\\ — —aH{a). Near the singular point, the op- 
tical depth can become large, thus invalidating the optically thin 
approximation. As a result, one cannot apply the popular equation 
(eq.|33ll to evaluate the brightness temperature, but should instead 
use the exact solution 

<5Tb(i.obs) = ar [Ts{^r){l - ^) - TcMB(ar)] [l " e-'-ob^] 

(47) 

where r^^.^^ is given by equation ( 145 1 when r^^i^^, > 1 or equa- 
tion ( 127b when t^^^^^ < 1 but not too small. The exact value of 
Ti/^jj^ is not important as long as t,^^^^^ ^ 1, since the r-dependent 
term should saturate 1 — exp {—r^^^^^ ) « 1 for large r^^^^ , in which 
case the brightness temperature is still finite, i.e. 



5Tb{iyobs) ~ ar [r,(e,)(l - -7) - rcMB(ar)] (48) 
instead of infinite as it would be using the popular equation i33\ . 



5.1.6 How Good is the Optically Thin Approximation during the 
EOR? 

It is often assumed that 21cm line is optically thin, fundamentally 
because 21cm hyperfine transition is highly forbidden with an ex- 
tremely small probability of Aw = 2.85 x lO^^^s^^. However, 
peculiar velocity gradients can enhance the optical depth in over- 
dense regions O To see this, we rewrite the optical depth in equa- 
tion l |27t as 



STt{zc, 



l + 5„ 



TCMB.O Q |1 + 5s,„| (1 - -^) 



= 0.00438 



0.02 

1 + (5pHI 



0.15 1 + Zcos ^HI,m('2cos) 



Qmh'^ 10 



0.5 



(49) 



where a = Ts (r) /Tcmb{zcos). For example in our reionization 
simulation at 2 --^ 9 when xm,m ~ 0.5, and with reasonable as- 
sumptions such as non-relativistic bulk motion v <^ c and small 
fluctuation 1 + Sp^^^ ~ 1, the optical depth can become of or- 
der unity when the velocity gradient is very negative, such that 
1 1 + Sq^v I < MM ^ as can happen in some overdense regions. This 
condition on the velocity gradient widens when 21cm occurs in ab- 
sorption (Ts < Tcmb) and becomes narrower when it occurs in 
emission (T^ > Tcmb). 

Figure [3] shows the PDF of the r^^^^ distribution of the 
IGM from our simulation data (see simulation details in § |6.3I >. 



" lUiev et alj i2003) showed that the 21cm line can become optically thick 
inside dense mini-halos, but this is a different effect from the enhancement 
due to velocity gradients we consider here. 



We smooth the N-body particle mass in the IGM onto a regular 
256'^ grid, compute the cell's velocity gradient using the SPH-like 
smoothing method described in Appendix [A] and compute the op- 
tical depth of the IGM. For simplicity, we drop the 1 — v/c factor 
in the optical depth calculation by assuming non-relativistic bulk 
motions. Figure[3ja) shows the extent by which velocity gradients 
alone can enhance the optical depth, by assuming Sp^^ = Sp^^. 
The PDF at large optical depth increases by roughly an order of 
magnitude when we decrease the spin temperature by an order of 
magnitude. In the pessimistic case (Ts/Tcmb = 0.1), as many as 
0.1% of the total cells have an optical depth of order unity. 

In Figure|3lb), we plot the same PDF for the actual Sp^^^ from 
the reionization simulation. Due to the "inside-out" character of 
reionization, the overdense regions that can have velocity gradients 
close to —1 ionize first and one can expect the effect to be much 
less. Figure[3fb) shows that for the case Ts/Tcmb = 0.1, only a 
fraction of up to 10^* of the total number cells approach an optical 
depth of 1. For the case Ts/Tcmb ~ 100 this fraction becomes 
as low as 10~^. We therefore conclude that we can safely use the 
optically thin approximation when calculating the 21cm radiation 
from the IGM. 

However, we should note that the optically thin approxima- 
tion may break down to a larger extent in one of the two following 
scenarios. 

(i) 21cm radiation from a halo or in spherical collapse at the 
turn-around point may be mostly optically thick, because Sg^y ~ 
— 1 there. The breakdown of the optically thin approximation may 
be more prominent when the 21cm line is in absorption against the 
CMB. 

(ii) When the 21cm radiation is computed directly from high- 
resolution particle data (and not from gridded data as above), a 
larger fraction of particles can be optically thick, since the parti- 
cle density is higher in overdense regions. 

The breakdown of the optically thin approximation merits more in- 
vestigation beyond the scope of this paper where we focus on 21cm 
radiation from the IGM and this approximation is mostly valid. We 
defer further analysis to future work. 



5.1.7 When the Mapping from Frequency to Position Along LOS 
is Multi-valued 

In the case of multiple 21cm transitions along the ray path, we la- 
bel the transition events by i = 1 , . . . , TV in sequence along the 
forward ray path. The optical depth starts from ro = (on the 
far side of the gas element), to Ti (after the ray passes through 
event i), and to rjv = T^ots- define the differential optical 
depth Ari = Ti — Ti^i — J . K^id^' which can be eval- 
uated using equation ( I27t with the transition location — > C»- 
To carry out the integration in equation (1181 . we split the inte- 
gral into a sum of A*' sub-integrals each over only one transition 
event, i.e., Jq""^" . ■ . dr^, — X^ili ^ ■ ■ ■ dT^, . Using the same 
trick as in § |5X3] we find Jo"°'"' 5'^,'(^')e"'^"'°>'.-<')dT^, = 
X^ili S^i {^i)ATi in the optically thin limit. Using the fact that 

'^''obB = J2iLi Ari, we find that (5T6(i^obs) = Y^iLi '^TbCfobs)!^,, 
i.e. the observed differential brightness temperature is the sum of 
contributions from all transitions, where each contribution can be 
evaluated using the equation for a single transition (eq|33t with r 
the position of each transition. 
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Figure 3. PDF of 21cm optical depth t,,^^^^ from our simulation data at z = 9.457 (50% ionized). The simulation is in a 114 Mpc/h box with IGM particle 
data smoothed onto a 256'' grid. PDF shows the probability of finding t^^j^^ in intervals of At^^j^^ = 0.1. We assume Tb/Tqmb = 0.1 (solid, black), 1 
(long-dashed, red), 10 (dotted, blue), and 100 (short-dashed, green). Left panel: assuming a fully neutral Universe (xjji = 1). Right panel: using the actual 
ionization pattern from the simulation. 



5.2 The Distortion of Apparent Location and Brightness 
Temperature by Peculiar Velocity 



In ij lS.ll we derived the equation for the observed 21cm brightness 
temperature, evaluating physical properties at the actual location 
r of the emitting neutral hydrogen atoms. However, observers can 
only determine the position of the source from the observed 21cm 
line frequency, i.e. in observer redshift space. To make theoretical 
predictions, it is therefore necessary to express the observed 21cm 
brightness temperature in observer redshift-space coordinates. This 
subsection deals with solving this issue. 

5.2.1 Distinguishing the Two Distortion Effects by Peculiar 
Velocity 

We should emphasize first that, although the effects of peculiar ve- 
locity on observed 21cm brightness temperature and on apparent 
location of sources are both due to the Doppler shift of the line 
frequency, the underlying mechanisms do differ. For the former, 
peculiar velocity distinguishes the emitter space from the FRW 
space, both of which are physical reference frames, and translates 
the difference between these two frames, through the transforma- 
tion of the line profile, to the optical depth that affects the bright- 
ness temperature measured by observers today in FRW space. It 
is a "real" effect in the sense that peculiar velocities change the 
observed brightness temperature, regardless of how observers in- 
terpret the location of source. 

For the latter, the observer redshift-space coordinates are sim- 
ply an artificial coordinate system that could be replaced by the 
observer real-space coordinates if observers could measure the pe- 
culiar velocities of sources and reconstruct the brightness temper- 
ature map in the sources' actual location. This is "artificial" in the 
sense that it is due to the observers' incomplete information on the 
location of the sources. 



The observed brightness temperature we derived in equa- 
tion l l36t is evaluated in terms of quantities measured in real space. 
We can rewrite equation ( 136) as 



1 



JCMB jar) 

r;''=«(r) 



(50) 



using the superscript "r" for real space explicitly. (See our conven- 
tion of superscripts "r" and "s" in § 13.11 ) By definition, the bright- 
ness temperature calculated from redshift-space quantities, ST^ (s), 
is equal to (r). So in principle, we can combine the two effects 
of peculiar velocity and find an expression for the brightness tem- 
perature using redshift-space quantities from 



(51) 



where r(s) is the inverse of the real-to-redshift-space mapping r 



g _ J. _l_ (i+Zobs) y f sjjow below that this relation can 
be simplified for the 21cm brightness temperature. We restrict the 
discussion to calculating the 21cm brightness temperature in the 
optically thin approximation here, since this is mostly valid in the 
IGM. 



5.2.2 21cm Brightness Temperature in Observer Redshift Space: 
Mathematical Approach 

We present the derivation of the equation for the 21cm brightness 
temperature in terms of redshift-space quantities in two ways. First 
in this subsection in a mathematical way, and in the next subsection 
in a more heuristic physical way. To simplify matters, we for now 
assume that Ts 3> Tcmb and generalize our result to an arbitrary 
Ts in fi \5JA\ 

Analogous to redshift space distortion in galaxy surveys, 
where the number of galaxies is preserved between real- to redshift- 
space, the number of emitting neutral hydrogen atoms is pre- 
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served in the 21cm signal, i.e. nfii(s)5V(s)(l + Zc 



) . From the real-to-redshift-space map- 
«ll (r) f , it is easy to find the relation 



n^ji(r)<51/'-(r)(l + 2cos 
ping, s = r + 'jri;^ 
between comoving volume elements in both frames 51/° (s) — 
5V''{r) |l + 5g^„(r)|. Therefore, the number density measured in 
redshift space is 



nm(s) 



(52) 



The mean number density must be preserved too, when averaged 
over a volume large enough to contain all gas of interest. In terms 



of fluctuations Sp (s) 



"Hl(s)-"Hl(Zcoa) 
nHl(2co3) 



where nni is the mean 



(physical) HI number density, we have 1 + Sp^^ (s) = | || , 

and hence in the Ts ^ Tcmb limit, 

Sn{s) = ^I^^nl,i{s) = ^tizcos) [1 + 5;^^{s)] . (53) 

This means that in the high Ts limit, the observed 21cm brightness 
temperature is directly proportional to the number density of neu- 
tral hydrogen atoms measured in observer redshift space. In other 
words, 21cm tomography maps exactly the neutral hydrogen distri- 
bution in redshift-space. This is the result we already found in §|4] 
but now more rigorously derived. 

In case of multiple transitions along the ray path, the bright- 
ness temperature is the sum of contributions from all transition 
events, as discussed in ij |5.1.7l Since these transitions correspond to 
the same observed frequency and therefore the same redshift-space 
location, equation l |53t still holds for the multi-transition case, since 
by definition the HI density in redshift-space is the linear addition 
of HI mass from all such transition spots per unit redshift-space 
volume. 



5.2.3 21cm Brightness Temperature in Observer Redshift Space: 
Physical Approach 

Now we rederive equation ( 153) by considering the physical mean- 
ing of brightness temperature. The 21cm brightness temperature 
is simply proportional to the specific intensity, i.e. 5Ti,{vo\^s) — 

2 

'2 — SIv where 57^ is the differential specific intensity rela- 
tive to CMB, and equal to the energy received from distant gas per 
unit observation time per unit transverse collection area per solid 
angle spanned by sources per unit observed frequency interval. The 
solid angle is proportional to the transverse area of the source, the 
observed frequency interval is proportional to the LOS distance in- 
terval in redshift space, and hence the energy received from a patch 
of sky near Vobs per unit time per unit collection area is propor- 
tional to the brightness temperature times the redshift-space vol- 
ume element. I.e., d^fi = dA''j_/d'A{ZohB), dUohB = |dS|| |/j/(Zobs), 

and dE/dtdAcou = C(zobs) '5^6(^'obs)<5V^^ where dA{zohs) is 
the comoving angular diameter distancj^ y{zohs) = Ao(l + 
Zobs)'^/-ff (zobs), dA"j^ is the comoving transverse area in red- 
shift space, cis|| is the comoving radial interval in redshift space. 



Here 



?0 



dA{z) = ^\n,\-y^S where E{z) 



-jj — is the relative cosmic expansion rate, and the function S(x) equals 



sin(a;) if Q,f^ < 0, x if Q.^, = 0, and sinhx if £7^. > 0. Strictly speaking, it 
should be (iA(^cos) that differs from dAC^obs) by (l + Zobs)/-'7(-2obs)- 
Since dA is large at high redshift, this difference is negligible. 



C{zabB) = 2fcsJ^obs/c^ciA(^obs)j/(zobs), and = dAljds|(| 
is the comoving redshift-space volume element. 

Consider a small region (e.g., a cell or a pixel) of 
the sky at the telescope's resolution scale. The detector 
simply smears subcell brightness temperature information 
by summing energies received from all unresolved sub- 
cells. For each subcell, STtSV 



5y(r)|l + 5S^„(r)| 
5Nm{r) = {1 + z. 



(1 + 



"HlCzcc 



l(r) 



where 



nui{f) 5V"^ (r) is the number of 
emitting neutral hydrogen atoms from the subcell at r. Ignoring 
the difference of observed frequency and redshift between the 
subcells, the brightness temperature of the cell is STtiiyohs) ~ 

^^(.(zcos)^^^ = ^^^(^co.) [1 + 5;hi(s)] (i-e- eq.m, where 
AV is the total redshift-space volume of the cell, and nfji ^.(,11 — 
(1 + zcosf {Y.SNm,suh) /AV = (1 + ZcosfANni/AV is 
the cell-wise (physical) HI number density in redshift space. 



5.2.4 Spin Temperature Reloaded 

In this subsection we generalize our calculation to the case 
of arbitrary spin temperature. Following the same algebra 
as in § 15.2.31 for each unresolved subcell, STt SV — 



temperature of a cell is 

5Tb{z, 



[1 - ^^^fIhI;^] • Then the brightness 



5Tb{uohs) 



nHl(Zcos) 



Tcmb 



T^cff 



(54) 



where 



nm 



Tcmb 



X 



1 



subcells 
TcMb(2cos 



E r^Hi(s) 



(55) 



sub 



is the redshift-space-volume-weighted cell-wise average of 

"-HI 

5Nj 



HI 



1^1 — '^y ™ j per unit proper redshift-space volume. Here we 
implicitly assume that spin temperature is preserved from real- to 
redshift-space, i.e., T^'^^is) = rj''=**(r). 



5.2.5 Breakdown of the Analogy to Galaxy Surveys 

From our results it is clear that the analogy to galaxy redshift sur- 
veys breaks down due to two effects: finite optical depth and finite 
spin temperature, as mentioned before in Sec.l4lF^ 

For the first case, when the IGM is optically thick to 21cm ra- 
diation, i.e., r^^ba ^ 1' brightness temperature is not linear in 



^® There is a third, more technical, difference between galaxy redshift sur- 
veys and 21cm surveys. In principle, the apparent location shift from real- 
to redshift-space results in the difference in the comoving transverse area 
and, hence, affects the redshift-space volume, in addition to the effect due 
to the change in the comoving LOS distance interval. This additional effect 
is non-negligible for galaxy redshift surveys at low redshifts, but small for 
high-redshift 21cm surveys (as discussed in Footnote |15t . We thank Antony 
Lewis (201 1, private communication) for pointing this out to us. 
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r^^jj^ (see eq. I47t . and the optical depth itself is affected by pecu- 
liar velocity through its dependence on spatial derivatives that are 
higher order than dv\\/dr\\ (see eqs.|45|and|46|l. Consequently, the 
brightness temperature is no longer proportional to the neutral atom 
density in redshift space. 

For the second case, e.g. at high redshifts where Ts 2> Tcmb 
is not satisfied^, neutral atoms in the same redshift-space volume 
element contribute unequally to the brightness temperature due to 
their spatial variation in level population, i.e., emitters can have dif- 
ferent luminosity. Thus the brightness temperature is no longer pro- 
portional only to the neutral atom density in redshift space. When 
the mapping from real- to redshift-space is single-valued, the pro- 
portionality between observed brightness temperature and neutral 
atom density in redshift-space is spoiled by the spatially-vaiying 
correction factor, 1 — TcMB/2T''^''^(r), according to equations ( I54t 
and l |55l l. However, in the more general case in which the mapping 
may be multi- valued, this coiTection factor is an average over the 
different real-space streams that contribute to the same redshift- 
space element, weighted by their different redshift-space neutral 
atom densities. 



5.3.1 Fully nonlinear power spectrum with finite optical depth 

Consider a slice 5T^(s) of a 3D data cube near Zcos, in red- 
shift space. The brightness temperature in Fourier redshift space is 
5T^{\i) = J rf^s e"* '' ST^{s). Since predictions of power spec- 
tra from theoretical modeling are made in real space, we should 
relate this to real-space quantities. The redshift- and real-space co- 
ordinates are related by equation (|6), and so the volume elements 
are related by d'^s — d^r 1 1 + Sg^^ (r) | . The observed brightness 
temperature is preserved (see eq. 1511 . and, in the general case of fi- 
nite optical depth, evaluated using equation J47b with optical depth 
using equation i49\ . The exact Fourier transform of brightness tem- 
perature in redshift-space is, 



5T^{k) = J d'r e-* "'- 
xTcMB.o 11 + <5s,,4r)l [a(r)( 



exp 
1 - 



4-1] [ 



I + Zc 
H{Zco 

II [1- 



fc||U||(r) 
-""Ob.! ,(56) 



where fcy = k ■ r. Note that t^^^^ is an implicit function of r, 
too. The fully-nonlinear power spectrum can be calculated by its 

definition /5T5f*(k)57f(k')) = (27r)^PAr(k)(5(3)(k - k'). 



5.3 Redshift-space Distortion on 21cm Power Spectrum 

The 21cm redshift-space-distorted p ower spectrum in the linear ap- 
proximation was explored in Barkana&LoebldlOOSi) . who showed 
that the linear 21cm power spectrum is distorted in a form anal- 
ogous to the linear redshift space distortion in galaxy surveys. 
The authors computed the power spectrum of linearized peculiar- 
velocity-corrected 21cm brightness temperature, nevertheless, in 
real space, i.e. they linearized gas density, neutral fraction, and 
particularly the velocity gradient correction 1/(1 + (5g^^(r)) ~ 
1 — (5g^„(r) by assuming 5q^^ <^ 1, and computed the Fourier 
transform of the brightness temperature evaluated in real space. The 
observable power spectrum, however, is in redshif t space. Although 
the ex pression of power spectrum derived in iBarkana & Loebl 
( I2OO5I) can give correct values on large scales, this approach is con- 
ceptually incomplete. In addition, the assumption of 5g^^ <C 1 
may break down on sm all scales. A further complication is that 
iBarkana & Loebl ( l2005h assume that the product of neutral frac- 
tion fluctuation and the gas density fluctuation, Sx^iSp^i, can be 
neglected, which can be invalid and cause the power spectrum to 
become inaccurate with a fractional error at the 200% level on small 
scales when the universe is 50% ionized ( Lidz et al. 2007) . 

In this section, we present a reformulation for computing the 
21cm power spectrum in observer redshift space, taking into ac- 
count both distortions in brightness temperature and in apparent lo- 
cation, and give the general equation for the linear redshift-space- 
distorted power spectrum without assuming either 5g^^ <^ 1 or 



It is generally assumed that sufficiently late after the fomiation of 
the first stars, the spin temperature is well above the CMB tempera- 
ture. This assumes, e.g., that the IGM is heated but only weakly ion- 
ized, as by the X-rays expected from early galaxies and mini-quasars (e.g. 
IChen & Miralda-Escudj liooi ). It also assumes that the first stars pro- 
duce a strong enough Lya-pumping background to couple Ts to the ki- 
netic temperature of the gas through the Wouthuysen-Field effect (e.g. 
ICiard i & Madau 2003). However, the length of the transition period from 
Ts < TcMB in the Dark Ages to Ts 3> T^CMB durin g the later stages of 
the EOR is an unsettled topic (see, e.g.. lBaek et alllOlQ) . 



5.3.2 Nonlinear power spectrum in the optically-thin 
approximation 

In the optically thin limit, we can use the approximation 1 — 
e^^^obs = Tvabs - before, the velocity gradient corrections for 
the optical depth and the redshift-space volume element cancel in 
equation l l56t . and we find that the fully nonlinear Fourier trans- 
form of brightness temperature in redshift-space in the optically 
thin limit is given by 

57f (k) = 5r,(z_) / d\ e'*- [1 + s;^^{r)] 



X exp 



1 + Zc 
H{Zco 



fc||«ll(r) 



Phi 

TpMB (^cos) 



(57) 



5.3.3 Quasi-linear fi]^-decomposition Scheme 

We work out a "quasi-linear" case in this subsection. In this we 
only take the density and velocity fluctuations to be linear, but the 
reionization fluctuations are allowed to be nonlinear. This means 
that we do not assur ne ^zh, <C 1 and thu s our approach is more 
general than that of IBarkana & Loe3 ( l2005h . We therefore choose 
not to call it "linear theory", but instead introduce the new name 
quasi-linear jiy^-decomposition scheme. 

On large scales corresponding to small enough k so that 
( ) '''II ^11 ^ ^' linearize the exponential and 

keep the linear term in v. We also linearize the spin- temperature- 
dependent term 



r,'-(r) 



TCMB (^cos 

rj'°«(r) 



(58) 



by defining its fluctuations as S^^ir) = [r?''(r) — f]{zcos)] /f]{zcos) 
where f7(zcos) is the mean value of 77. We keep only the linear terms 
in velocity, neutral density fluctuations, and T^-fluctuations, and 

* ''a''(r) is the Fourier 



find <5T,'''''""(k) = STt{z,o.)f]{zco.) 



Here a'- (k) = J d^ 



transform of the quantity a'^(r) in real space. On large scales, the 



© 201 1 RAS, MNRAS Q00.[Tlt30l 



21cm Redshift Space Distortion: Methodology Re-examined 1 5 



velocity field is linear, U||(k) = i y ]^^^°°° j Sp^(k.)^, where 

= fc = |k|, and Sp-^Ck) is the total hydrogen density 

fluctuation in Fourier real-space. So we find 

(59) 

The power spectrum in the quasi-linear /ik -decomposition scheme 
in redshift space, defined as ^^T^'""" {k)ST^-'^^''\k')^ = 
(27r)='P^'^""(k)5(=''(k-k'),is 

Pl'^""(k) = P^„ (fc) + Pp2 {k)^i + P^4 (fe)Mk , (60) 

where the moments of /ik -polynomial expansion are 

+2n;„,,.,(fc)] , (61) 

p,2 (fc) = 2 (5r,77) ' [p,;^^ (fc) + p,; (fc)] , (62) 
p^4(fc) = (5r,77)'p,;^,,^^(fc), (63) 

where all quantities here depend implicitly on the redshift Zcos- 
Here P^^ denotes the auto power spectrum of the quantity a^{r), 
and P^ f, is the cross power spectrum between fields a' (r) and 
b^{r), both in real space. Note that, strictly speaking, the power 
spectra involving 5^ are not statistically isotropic due to the distor- 
tion by peculiar velocity as in equation ( I34t . Since the correction 
is of order v/c, we ignore it here. When Ts 3> Pcmb , rj — 1 and 
Srj = 0, and the power spectrum in quasi-linear /ik -decomposition 
scheme reduces to equation (|4}. 

Although we derived the scheme by assuming linear density 
and velocity fluctuations, when using it on simulation data, we nor- 
mally use the non-linear density fluctuations given by the simula- 
tion. 

As pointed out above, each moment of the /ik -decomposition 
can contain higher-order auto- and cross-correlations involving 
density and ionization fluctuations, because Sp^^ = 5p^ + Sl-^^ + 
^PH ^xai- To see this explicitly, for example, in the simple case 
Ts ^ TcMB in which r} = 1 and 5,, = 0, we can rewrite the 
moments as follows. 

p^2(fc) = 2 5r'p,;^^.,^^(fe) 

= 2 5T^P,;^,,^^(fc)+P,:^^,,^^(fc) 

+-P^^hi*ph.*.h('=)] ' (65) 

P^4(fc) = STlPl^,,^^{k) . (66) 

However, the quasi-linear /^k -decomposition scheme neglects 
the nonlinear coupling of peculiar velocity and io nization fluc- 
tuatio ns, which we will investigate in future work JShapiro et al.l 
l20Tlh . 



5.3.4 Linear Theory 

lBarkana&Loebl ( l2005h linearizes both density and ionization fluc- 
tuations, and discards all three- and four-point correlations in the 
expansion of moments, i.e. in the simple case Ts S> Tcmb, equa- 
tions ( 164166b reduce to 

p^o(fc) = 5r^p,:^^,,^^^(fc) + 2Pi:^^,,^^(fc) 

P^. (fc) = 2 5fl [P,;^ (fc) + Pl^^ (fc)] , (68) 

p^4(fc) = 5r'p,;^,,^^(fc). (69) 

iLidz et al. demonstrated that, if peculiar velocity is not 

taken into account, i.e. only zeroth moment is concerned, the ne- 
glect of higher-order correlations can result in significant errors in 
21cm power spectrum. They also pointed out that, for the same 
reason, 21cm redshift-space power spectrum computed using the 
linear theory of Barkana & Loeb (2005) can have large errors, but 
they did not provide any detail or analysis of computing the non- 
linear power spectrum, nor did they propose an analytic solution 
that incorporates all of the relevant higher order terms. 

In our paper, in addition to investigating the fully nonlinear 
power spectrum, we propose the quasi-linear /^k -decomposition 
scheme as a solution that can as well separate the cosmological 
density fluctuati ons from the ioniz ation fluctuations just as the lin- 
ear theory (Bark ana & LoeblboOSh does, but account for higher or- 
der correlations due to nonlinear ionization fluctuations. 

6 COMPUTATIONAL SCHEMES TO PREDICT 

BRIGHTNESS TEMPERATURE IN REDSHIFT SPACE 

6.1 Exact Steps in the Case of Finite Optical Depth 

Analytical models and semi-numerical or numerical simulations 
provide us with rea/-space data. In order to make predictions for 
the observed 21cm power spectrum, we need to calculate the fully 
nonlinear 21cm brightness temperature accurately and efficiently 
in redshift space, accounting for all effects of peculiar velocities. 

As explained in § 15.2.11 the effects of peculiar velocity can 
be separated into an effect on the observed brightness temperature 
and one on the apparent location of the 2Icm emission source. So 
in principle, in order to compute the signal in redshift space, the 
brightness temperature should (1) first be corrected by the velocity 
gradient, evaluated in real space, using the exact formula of 21cm 
brightness temperature (eq. 147 ) with finite optical depth (eq. |49l l, 
and (2) then shifted to the apparent location corresponding to the 
Doppler frequency shift, with the volume element re-sized accord- 
ing to the velocity gradient, and (3) finally resampled onto a regular 
grid in redshift space. Power spectra calculated this way should be 
equivalent to those using equation J56t . This process is in general 
computationally cumbersome. 

6.2 Real-to-Redshift-Space-Mapping (RRM) Schemes 

Since the optically-thick cells are very rare in the IGM, as we 
have shown in § 15.1.61 we may evaluate brightness temperature in 
the optically-thin approximation (eq. [36}. In doing this, although 
brightness temperature in an optically-thick cell would become ar- 
tificially divergent in real-space, its net contribution to the bright- 
ness temperature in redshift-space is still finite and proportional to 
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the total number of neutral atoms in that cell, because the redshift- 
space volume element of this cell is compressed accordingly. This 
has been well discussed in § 15.21 We can exploit the proportional- 
ity between the 21cm brightness temperature and the neutral atom 
number density both measured in redshift space. Inspired by com- 
mon wisdom in large-scale structure simulations, we propose two 
computational schemes based on mapping the neutral atom density 
from real- to redshift-space, and then computing the 21cm bright- 
ness temperature in redshift space using equation ( 153V We also as- 
sume Ts ^ TcMB in this section, but our schemes can be readily 
generalized to the arbitrary Ta case. 

Strictly speaking, these two schemes are accurate only when 
the optically-thick cells are rare enough, because neutral atoms in 
those cells should be "self-shielded" to 21cm radiation. We will 
revisit in detail the accuracy of power spectrum in the optically- 
thin approximation in § |7] 

6.2.1 Partide-to-Partide-to-Mesh(PPM)-RRM scheme 

Most numerical simulations of reionization are processed as fol- 
lows. First one runs a large-scale N-body simulation, from which 
one obtains gridded density fields and the collapsed halo informa- 
tion such as location and mass. The reionization simulation is then 
run on these gridded density fields using the halos as sources of 
radiation. Since the RT grid resolution is typically coarser than the 
N-body particle resolution, the most accurate 3D map of the neu- 
tral atom distribution in redshift space that can be possibly achieved 
from a given reionization simulation is made by taking advantage 
of the high-resolution N-body particle information. We propose 
the Particle-to-Particle-to-Mesh Real-to-Redshift-Space-Mapping 
("PPM-RRM") scheme as follows: 

• We compute the bulk-flow velocity of the IGM at the position 
of particles directly from N-body particle data using an adaptive- 
kernel, SPH-like approach. The SPH-smoothed bulk velocity as- 
signed to each particle is the smoothed momentum density divided 
by the smoothed mass density, evaluated at the particle location. [3 

• We assign each particle the neutral fraction from the RT grid 
cell that it is located in. 

• For a given LOS direction, we Doppler-shift the N-body par- 
ticles to their apparent locations according to the LOS bulk-flow 
velocity, in accordance to equation (|6}. 

• We compute new smoothing kernel lengths using the new par- 
ticle positions in redshift-space. 

• We use those kernel lengths to smooth the particle data (i.e. 
H I mass) onto a regular, redshift-space grid (see the discussion of 
grid resolution below). In this step, we exclude particles contained 
in haloQ 

• From this latter, gridded density field, we compute the H I den- 
sity fluctuations in redshift-space, and from this the 21cm bright- 
ness temperature measured in redshift space using equation i53\ . 

If a hydrodynamical simulation is coupled to N-body cold dark matter 
(CDM) simulation, then the gas particle velocity can be directly used. But 
since our simulations are dark matter only, we approximate the gas bulk- 
flow velocity as the SPH-smoothed velocity at the particle location (see 
Appendix IaI. One cannot use the particle velocities directly because those 
can be multi-streaming. In all this we assume that the gas traces the dark 
matter exactly, which is a good approximation on large scales. 

We simulate the reionization of the IGM, and therefore compute the 
21cm brightness temperature only from the IGM, so excluding particles in 
halos. 



Some details of the particle smoothing algorithm are discussed in 
Appendix IaI We use adaptive kernels rather than fixed-kernels so 
as to better resolve the small scale spatial variations in overdense 
regions. 

The high wavenumber modes in the power spectrum can be in- 
accurate due to sampling effects when calculating the power spec- 
trum using the fast Fourier transform (FFT). Instead of correcting 
the power spectrum using the method proposed by Jing, (2005]), 
we partly avoid the sampling effect by gridding the particle data 
onto a redshift-space grid at four times higher resolution than the 
RT grid, but only keeping the modes in the power spectrum with 
k ^ n/AL (AL is the RT grid spacing), i.e. one-quarter of the 
Nyquist wavenumber for a grid with the resolution AL/4. The rea- 
son for this and a summary of the sampling effect are discussed in 
more detail in § 16.41 

The PPM-RRM prescription can be summarized as "Pr — > 
Ps — )■ Ms(4x RT)" where "P" means particle data, "M" means mesh 
data, subscript "r" means real-space, "s" means redshift-space, and 
"4xRT" indicates that the grid resolution is 4 times finer than RT 
grid resolution. Figure |4] shows the flow chart for the PPM-RRM 
scheme. 



6.2.2 Mesh-to-Mesh(MM)-RRM scheme 

Manipulating N-body particle data is accurate but computationally 
costly (see Table [T] below). Since the N-body particle data typi- 
cally already have been smoothed onto a regular, real-space grid 
in order to simulate the radiative transfer, we propose an alterna- 
tive scheme, the Mesh- to-Mesh Real-to-Redsh ift-Space-Mapping 
("MM-RRM") scheme. iMellema et all t2006hl ) were actually the 
first to use the MM-RRM scheme to produce brightness tempera- 
ture spectra and maps along the LOS (their Figs. 4, 9 and 10), but 
did not provide a detailed description of the method in their paper. 
This scheme saves computational resource by using the real- space 
grid data such as cell-wise mass density, velocity, and ionization 
fraction, but gives consistent results (depending on the grid resolu- 
tion, to be tested in §|6j6j. The MM-RRM scheme works as follows: 

• As the preliminary step, we grid the N-body particle data in 
the IGM (i.e. particles in the halo excluded) onto a regular, real- 
space grid with a resolution n times finer than the RT resolution, 
using our adaptive kernel SPH-like smoothing. This provides us 
with cell-wise density and velocity fields. 

• We assign each cell the neutral fraction from the RT grid that 
this fine cell belongs to. 

• We assume the cell-wise velocity to be the velocity at the cell 
center, and compute the LOS velocity at the boundary between two 
LOS-neighboring cells by linear interpolation. 

• We shift the cell boundaries to their apparent locations ac- 
cording to their LOS velocity, in accordance with equation (|6]l, 
whereby the real-space cell can get stretched or compressed in red- 
shift space. In high density cells the boundaries of a cell can cross 
each other in redshift space, an effect known as the finger of God. 
When this happens, we switch the cell's crossing boundaries so that 
the cell size is always positive. 

• We regrid the neutral hydrogen mass from the real-space grid 
onto a regular, redshift-space grid at the same resolution, by count- 
ing the overlapping volumes; e.g., if the LOS is along the a;-axis, a 
real-space cell k) with the size Ax stretches to the length Ax[ 
in redshift-space, with a portion of this length, AL^y , overlap- 
ping the cell (i'j j, k) in the regridded, redshift-space, mesh, then 
all real-space cells (i, j, k) contribute to the neutral hydrogen den- 



© 201 1 RAS, MNRAS 000.[Tlt30l 



21cm Redshift Space Distortion: Methodology Re-examined 1 7 



Preliminary Data Processing 



N-body Particle data (x,v) 



RT Grid Data (xj) 



Compute Bulk-flow Velocity 



Paint Neutral Fraction on Particle Mass 



Particle-to-Particle-to-Mesh Real-to-l'edshift-Space IVIapping 



Set the Line of Sight (LOS) Direction 



Doppler Shift Particles to Redshift-space Locations by LOS Bulk-flow Velocity (Particle-to-Particle) I 



Find New Smoothing Kernel Lengths 



Resmooth Halo-excluded Neutral Hydrogen Mass onto a Regular Redshift-space Grid (Particle-to-Mesh) 



Figure 4. Flowchart of the PPM-RRM scheme. 



sity of the redshift-space cell {i' ,j, k), according to 

Pm{i',hk) = ^Fi.i, p'iii{i,j,k) , (70) 

i 

where j/ is the fractional overlap of the real-space volume i with 
the redshift-space volume i' , i.e. F^ ^/ = ^i^i.v I (the indices 
J and k are not relevant here because we move all cells along the 
X-axis). 

• We compute the HI density fluctuations in redshift-space, and 
from this the 21cm brightness temperature using equation (153b . 
This is done at at n times higher resolution than the RT grid, but 
when calculating the power spectrum we only keep modes with 
k ^ Ti / /\L (AL is the RT grid spacing). 

The MM-RRM scheme can be summarized as "[Pr 
Mr(nxRT)]— >■ Ms(nxRT)", where the operation inside the square 
bracket is the prerequisite step. In ij |6.6| we will experiment with dif- 
ferent resolution factors n to find the optimal resolution. Figure |5] 
shows the flow chart for the MM-RRM scheme. 



6.2.3 The Redshift-space-distorted Lightcone Ejfect 



Both the PPM-RRM and MM-RRM schemes deal with simula- 
tion data from a finite volume at a fixed cosmic time, implicitly 
assuming that the cosmic evolution of both neutral fractions and 
density perturbations are negligible during the light travel time 
across the simulation box, f cross- For the typical simulation vol- 
ume sizes (100-200 Mpc) one does not expect much evolution in 
the density field during tcross ■ However, the neutral fractions may 
evolve much more rapidly during some periods of the EoR. If 
{d\nxi/dt)5t > 1, then we must take into account this so-called 
lightcone effect ('Bar kana & Loebll200 6^ and couple it to peculiar 
velocity. This implies first time-interpolating the particle data to 
the appropriate look back time and the corresponding real-space 
location and then shifting the particles to their apparent location ac- 
cording to its interpolated LOS peculiar velocity, and finally map- 
ping these time interpolated particles onto a regular redshift-space 
grid on the lightcone. The full version of the lightcone PPM-RRM 



© 201 1 RAS, MNRAS 000.[Tll30l 



18 Y.Maoetal. 
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ICM N-body Particle Data (x,v) 



f Coarsened ICM N-body Crid Data {cellwise x, v) J 



RT Crid Data (x_i) 



Coase Crid Real-space Neutral Hydrogen Ceilwise Density and Veiocity of ICM 



Mesh-to-Mesh Real-to-Redshift-Space IVIapping 



Set the Line of Sight (LOS) Direction 



Interpolate the Neighboring Cellwise Velocities to be the Cell Boundary Velocity 



Move Cells by Doppler Shifting Cell Boundaries to Redshift-space Locations (Mesh-to-Mesh) 




Resmooth Neutral Hydrogen Mass onto a Regular Redshift-space Crid by Counting Overlapping Volumes 



Figure 5. Flowchart of the MM-RRM scheme. 



scheme is beyond the scope of this papei[__| and we postpone an 
investigation of this effect to a future paper in this series. 

6.3 Simulations 

For our reionization simulation we use a new large-scale, high- 
resolution N-body simul ation of the ACD M universe (performed 
with the CubeP^M code, illiev et"ani2008bh in a comoving volume 
of Lbox = 114Mpc//i on each side using 3072^ (29 billion) parti- 
cles. To find the halos, we use the spherical overdensity method and 
require them to consist of at least 20 N-body particles; this implies 
a minimum halo mass of 10* Mq. 

Assuming that the gas traces the CDM particles exactly, we 
grid the density on a 256^ grid using SPH-like smoothing with an 
adaptive kernel. The halo lists and density fields are then processe d 
with the radiative transfer code C^Ray l lMellema et"an 12006 j) . 
Each halo releases f~, ionizing photons per baryon per At = 



"^^ 'Melle ma et alj j2Q06 bh did apply such a time interpolation of grid data, 
both on the neutral fraction and density fields. 



11.5 Myrs, with f-, = 150 (/^ = 10) for halos below 10^ 
(above 10^ Mq ), respectively. To incorporate feedback from reion- 
ization, halos less massive than 10^ A/q located in ionized regions 
are not producing any photons. 

The simulations were run on the University of Texas Sun Con- 
stellation Linux Cluster Ranger, one of the largest computational 
resources in the world. Both codes are massively parallel, using 512 
compute nodes, e ach with one Quad-Co re 64 -bit proces s or. We re- 
fer the readers to iFriedrich etall (l201ll) and Illiev et al.l (l201lh for 
more details of this simulation which in those papers is labelled as 
"163Mpc.g8.7_130S". 

The simulations used the following set of cosmological pa- 
rameters Ha = 0.73, = 0.27, = 0.044, /i = 0.7, erg = 
0.8, — 0.96 where Hq = 100/t km Mpc~ ^ , consistent with 
the WMAP seven-year results dKomatsu et alj2oTll) . 



6.4 Sampling Effects 

Measuring power spectra using a fast Fourier transform (FFT) of 
gridded data suffers from the so-called sampling effect. This effect 
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Figure 6. Aliasing elfect in the PPM-RRM scheme: 21cm redshift-space 
ID power spectrum at z = 9.457 (50% ionized), when the particle data 
is smoothed onto a regular, redshift-space, grid with the RT grid resolution 
(256^, long-dashed, blue), or four times finer (1024'^, solid, black). The 
vertical lines are at fe = fc^^^'/4 = 1.75h/Mpc (thin long-dashed) and 

k = fc^°^''V4 = fc*^^*^' = 7h/Mpc (thick dot-long-dashed), respec- 
tively. The fractional error plotted in the inset is with respect to the power 
from the 1024^ grid. 



is due to the mass assignment of particle data or continuous fields 
to a chosen grid. In cosmology, it was first extensively discussed for 
power spectru m measurements of density fields in large scale struc- 
ture (see, e.g., |jing||2005L l Cui et al.ll2008l and references therein). 
The mass assignment is equivalent to convolving the true density 
field with a window function and sampling this convolved density 
field with a finite number of grid points. The power spectrum of the 
convolved field is a biased one, i.e., (Jing,2005) 

P-^(k) = ^|w'(k + 2fcivn)|'p(k + 2fcivn) + P,hot, (71) 

n 

where (k) and P(k) are power spectra of the convolved and 
true field, respectively, W(k) is the Fourier transform of the win- 
dow function, Pgiiot is the shot noise, and the summation is over all 
three-dimensional integer vectors n. The sampling effects include 
three aspe cts that can affect the true power spectrum measurement 
jCui et al . 2008). 

• Smoothing effect: the Fourier window function |iy (k)p falls 
off shaiply from \W{0)\'^ = 1, e.g., for a Cloud-ln-Cell (CIC) 
window function, \W\^ = 0.90 at fc = fc]v/4, but \W\'^ = 0.66 at 
k = A:jv/2, where /cat = 7r/a is the Nyquist wavenumber for some 
grid spacing a. 

• Anisotropy effect: the Fourier window function is not isotropic 
for a given k, and the anisotropy is significant for ~ fc]v. 

• Aliasing effect: higher wavenumber modes (n 7^ 0) contam- 
inate the true mode at k, preventing us from relating P^(k) and 
P(k) straightforwardly. For a FFT, (— fc]v,fcjv) is the range in 
k-space that a finite resolution grid can probe. Thus those high- 



wavenumber modes that contaminate are due to modes of the unre- 
solved field below the grid resolution. 

The smoothing effect and anisotropy effect can easily be cor- 
rected for, e.g. by just deconvolving P^ (k) with the normaliza- 
tion |W^(k)p. Correcting the aliasing effect is more difficult, and 
may be done using the iterative m ethod proposed and tested for 
the density power spectrum by Jin3 J2005h . Instead, we can be less 
ambitious and define a "comfort" zone (k ^ some critical value) 
where the FFT power spectrum has negligible errors. This can be 
done because at low enough k, all these samp l ing ef fects should be 
insignificant. The test problems in both Jiijgl ( |2005|) and lCui et al.l 
l [200S) seem to agree that the raw density power spectra for differ- 
ent window functions agree at fc < fe]v/4. Here we test this on the 
21cm power spectrum. In Figure |6] we compare two power spec- 
tra, both calculated with the PPM-RRM scheme but differing in 
the resolution chosen for gridding the redshift-space particle data, 
256'^ and 1024'^, respectively. As can be seen in the figure, both 
power spectra agree for k < fe^^^' /4, where k^^'^^^ is the Nyquist 
wavenumber of the 256^ grid. We therefore conclude that if we use 
the 256'' grid, we can trust the results for k ^ fcj^^^' /4. 

However, this comparison also shows that we can use our 
high-resolution N-body data to try to capture the modes between 
k = fc^''*'' /4 and k = fc^''*'' . By sampling the Doppler-redshifted 
particle data onto a grid with a resolution of 4xRT = 1024'* we can 
minimize the smoothing and anisotropy effects. We also minimize 
the aliasing effect due to the gridded density and velocity data. The 
aliasing effect due to the finite resolution of the ionization frac- 
tion field obviously cannot be coiTected for this way. However, this 
effect may be quite small due to the nature of the ionization frac- 
tion field. Recall that the aliasing effect is due to the contamination 
from high-wavenumber modes unresolved by the grid resolution. 
For blackbody type sources, the edges of ionized regions are shaip, 
i.e. the ionization fraction is very close to 1 inside and very nearly 
outside ionized regions. Therefore only cells that contain bound- 
aries of ionized regions have unresolved subcell information. The 
fraction of boundary cells for an ionized region of A'^ cells in each 
dimension is ~ N^/N^ = l/N. The peak of the H II bubble size 
distribution can be ~ 10 Mpc, corresponding to ^ 22 RT cells 
across a bubble (see, e.g. iFriedrich et al.ll20lll) . For such bubbles 
only ~ 4% of the cells contribute to the aliasing effect. Only if 
there are many small bubbles of size less than an RT cell, would 
the ionization field introduce a substantial aliasing effect. 

Given this argument it would seem prudent to choose the 
smoothed grid resolution to be four times smaller than the RT res- 
olution, as this minimizes the sampling effects for the 21cm power 
spectra. We therefore adopt this approach. The modes between 
fcjv/4 and fcjv (where fcjv here corresponds to RT grid resolution) 
may still be affected by the aliasing effect due to the finite RT grid 
resolution, but we expect this to be a minor effect. 



6.5 Tests of PPM-RRM Scheme 

Since the PPM-RRM (4xRT) scheme retains the particle data the 
longest by mapping them directly into redshift space, it can be ex- 
pected to be more accurate than the MM-RRM scheme. We there- 
fore first present tests for the PPM-RRM scheme in this section and 
in the next section compare the results of the two schemes. 
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6.5.1 Conservation of Mass 

The mean total (and neutral) hydrogen density is conserved be- 
tween real- and redshift-space, because (i) the total (and neutral) 
hydrogen atom number is conserved, and (ii) the the total space is 
conserved for a volume large enough (peculiar velocity vanishes 
for large distances) or a periodic box, since J SV^Sg^^ (r) is a total 
derivative. 

For the simulation box, the total volume is automatically con- 
served. We can therefore check whether our schemes conserve 
mass by checking the conservation of mean hydrogen density and 
HI density. Conservation of the mean density could be violated 
if a scheme would undercount particles after shifting particles to 
redshift-space. 

Our code passes this test by showing that the fractional dif- 
ferences of the mean total (and neutral) hydrogen density between 
real-space and redshift-spaces with three distinct LOS directions 
are zero, i.e. smaller than machine error. 



6.5.2 Large Scale Test 

As shown in § 15.31 the fully nonlinear power spectrum reduces to 
the quasi-linear /ik-decomposition power spectrum at large scales. 
We use this here to test the PPM-RRM scheme. Figure|7](top pan- 
els) shows the ID dimensionless power spectrum A2icm(^) = 
k'^ P^j^(k) /2tt^ calculated with the PPM-RRM scheme. In order 
to minimize noise, we averaged the power spectra for three distinct 
LOS directions (namely, along x-, y- and z-axes). Plotted in the 
same figure is the ID quasi-linear ^k-decomposition power spec- 
trum calculated directly from the real-space ionization fraction (on 
the 256'^ grid) and density and velocity data (on the 1024'^ grid), 
using equations l |60t - i63\ . We choose the 50% ionized epoch for 
this comparison. Note that even though we use the quasi-linear /ik- 
decomposition scheme equations to evaluate the power spectrum, 
we use the fully nonlinear density and ionization fraction fields 
from the simulation. 

The comparison shows that the nonlinear power spectrum 
computed from the PPM-RRM scheme agrees with the quasi-linear 
^k-decomposition power spectrum at large scales (k < 0.3 h/Mpc) 
within 5%. This confirms that the 21cm brightness temperature data 
cube constructed by the PPM-RRM scheme captures the correct 
large-scale fluctuations in redshift space as dictated by the quasi- 
linear /ik-decomposition scheme. The nonlinear power spectrum 
deviates from the quasi-linear /ik-decomposition power spectrum 
at intermediate scales (0.3 < fc < 2 /i/Mpc) at the level of ~ 10%, 
and even larger deviations can be found at smaller scales. In the 



To parallel-process N-body particle data using Message Passing Inter- 
face (MPI) software, the simulation volume is divided into cubic partitions 
and particles are assigned to the partition within which they are located. 
Each partition is processed independently by a given node in the parallel 
computer. The mapping described here of particle locations from real- to 
redshift-space can move a particle out of its original (real-space) partition 
into another, even to one which is not a neighbor partition. In that case, the 
number of partition pairs that must share particle data, to exchange parti- 
cles, can be large and, hence, computationally inefficient. Fortunately, we 
find that the size of each partition in our N-body simulations (which we 
also use for our real- to redshift-space mapping) is larger than the maxi- 
mum Doppler shift of particles in comoving coordinates, so only neighbor- 
ing partitions need exchange particle data. 

It still has the unit mK^ . It is dimensionless with regard to Fourier space 
units. 



second paper of this series ( IShapiro et al.ll20lTl) . we will investi- 
gate in detail the cause of this departure from linearity and how it 
affects the use of 21cm observations for cosmology. 

6.5.3 Test Down to Small Scales 

Similar to the quasi-linear /ik-decomposition scheme test that em- 
ploys the real-space grid data to compute the redshift-space statis- 
tics, we can compute the redshift-space power spectrum at all 
scales, in principle, by evaluating the integral in equation J57I ). The 
integration can be carried out by a fast Fourier transform of the 
data cube F(r) = exp [-i (^i^) fc||«|| (r)] • [l + S;^^{r)] 
(assuming Ts 3> Tcmb) for any given /cy, and then picking up 
only those modes with the LOS component i.e., 5T^(k) = 

STb{zcos) F(k) only if k ■ n — fcy. We can construct the whole 
Fourier data cube by making such FFT evaluation for each k-space 
plane of constant fc|| ^ 0, exploiting the symmetry 5T^{—]i) — 

ST^ (k), with fc[| discretized in units of 27r/Lbox. 

Note that in order for the discrete Fourier transform to 
be a good approximation to the continuous Fourier trans- 
form, the particle data should in principle be smoothed to 
compute the cell-wise average {F{r))^^^^ and in particular 

<^exp ^—i ^ hIz""") ) ^^ll^ll ^ directly. However, to take ad- 
vantage of existing cell-wise density and velocity data on 
the grid, we evaluate ^exp i ^i^^^^ ^H'^'ll ^ — > 

'^^P ("StStj) ^11 (^ll('"))ccii]- compute the ID power 
spectrum from the Fourier modes, averaged over three independent 
LOS directions. We name this method of evaluating the power spec- 
trum the Direct Evaluation by Multiple Real-space FFTs (DEMRF) 
scheme. 

Obviously, the DEMRF scheme is accurate only when the cell 
size is not too small so that (v'^iv)^ ~ (v^\ ('^))ccii ^ 

1 in the Taylor expansion of ^exp i (^-^y?^) ^II'^II ) 
On the other hand, if the grid is too coarse, the high-fc powers are 
subject to the sampling effect and become inaccurate. We experi- 
ment on the trade-off by trying out the DEMRF scheme on grid data 
with different resolutions (256^, 512^ and 1024^), and find that for 
a box of size 114 Mpc//i, the 1024'' grid is too fine and fails to 
make sensible results due to the subcell nonlinearity. We plot the 
DEMRF result computed from 256"' and 512'' grids in Figure |7] 
(bottom panels), and find that within the comfort zone for each grid 
(1.75 Mpc//i and 3.5 Mpc//i, respectively), the PPM-RRM result 
agrees with the DEMRF results within 1%. 

The three tests presented thus show that the PPM-RRM 
(4xRT) scheme is accurate on both large and small scales. We can 
now use this to test our other scheme. 

6.6 Test of MM-RRM Scheme 

The MM-RRM scheme is expected to be less accurate than the 
PPM-RRM scheme since it grids the particle data before moving 
to redshift space and inevitably small scale information is lost in 
the process. For example, the gas density within an RT cell is as- 
sumed to be uniform, so that the resized cell in redshift-space can 
be uniformly regridded by counting overlapping volumes. This as- 
sumption obviously ignores the subcell dumpiness. Second, the 
scheme treats the velocity of cell boundary as the linear interpo- 
lation between cell-wise velocities of two neighboring cells and 
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Figure 7. Tests of the PPM-RRM scheme for the 21cm redshift-space ID power spectrum at 2 = 9.457 (50% ionized). Top panels: large scale test of the 
PPM-RRM scheme (solid, black) against the quasi-linear ^^-decomposition scheme (long-dashed, red), both computed from a 1024^ grid, and Fourier modes 
kept only at k ^ A:^''^^-'/4 = fc^""^' = 7h/Mpc (the thick vertical lines). The fractional error is that of the quasi-linear /i^ -decomposition scheme result 
with respect to the PPM-RRM result. 

Bottom panels: small scale test of the PPM-RRM scheme computed from 1024'' grid (solid, black) against the DEMRF scheme, computed on a 256'^ grid 
(dotted, green) and on a 512^ grid (dot-short-dashed, blue), respectively. We only keep modes k ^ ^j^''^' = 7 h/Mpc. The black and blue curves are 
almost indistinguishable until at the very large k in the left bottom panel. The fractional errors are with respect to the PPM-RRM result. The vertical lines at 
k = fc'^^'^V^ = 1.75 h/Mpc (thin long-dashed) and k = fcj^^^'/4 = 3.5 h/Mpc (thick dot-long-dashed) delimit the comfort zone for the D£M/?F result 
computed on a 256^ and 512'^ grid, respectively. 



thus ignores small-scale velocity fluctuations at the inter-cell scale. 
Third, the treatment of cell boundary crossing is approximate and 
a careful treatment should require particle data to mimic the finger 
of God effect. However, the scales that are affected depend on the 
resolution chosen, and if one can choose a grid with fine enough 
resolution, there is a hope that the MM-RRM scheme can yield as 
accurate power spectra at fc ^ fcjv (corresponding to RT grid reso- 
lution) as the PPM-RRM scheme does. 

We experiment with the resolution of the MM-RRM scheme 
by choosing n = 1, 2 or 4 in the pipeline "[Pr — ^ Mr(nxRT)]— > 
Ms(nxRT)" (where nxRT means the grid resolution n times finer 
than RT grid resolution). We compute the 21cm power spectrum 
for each of these three resolutions and plot them for the modes 



k ^ in Figure [8] As above we average over three LOS 

directions. We use the PPM-RRM(4xRT) result as a benchmark. 
All MM-RRM results agree with the PPM-RRM result down to 
the scale fc < 1 h/Mpc, while at high k the MM-RRM(1 xRT) and 
(2 X RT) results deviate from the benchmark by up to 40% and 20%, 
respectively. Fortunately, the MM-RRM(4xRT) result agrees with 
the benchmark within 1% error on all scales down to k ^ fc^^^'- 
We therefore conclude that MM-RRM(4xRT) gives as accurate re- 
sults as the PPM-RRM(4xRT) scheme for fc < fcjy. 
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Table 1. Usability, accuracy and efficiency of vaiious computational schemes for the redshift-space brightness temperature. Our simulation is in a box with 
1 14 Mpc//i on each side, has 3072^ N-body particles, and evolves reionization on a 256^ RT grid. 



PPM-RRM MM-RRM DEMRF Quasi-linear 

fH^ -decomposition 

TxRT 4xRT TxRT 2xRT 4xRT TxRT 2XRT 4XRT 



Input data type 



N-body particle (x, v) 
(3072^ particles), and RT 
grid Xi (256^ grid size) 



Cell- wise (x, v) in 256^, 
512^, and 1024^ grid size 
(Ix, 2x and 4xRT, respec- 
tively), and RT grid Xi (256^ 
grid size) 



Cell-wise (x, v) in 256^ 
and 512'^ giids (Ix and 
2x, respectively), and RT 
grid Xi (256^ grid) 



Cell-wise (x, v) in 
1024'^ grid size, and 
RT grid Xi (256^ grid 
size), or real-space 
power spectra directly 



Output data type 



HI density in redshift- 
space grid in 256^ and 
1024^ grid size (Ix and 
4 X RT, respectively) 



HI density in redshift-space 
grid in 256^, 512^, and 1024^ 
grid size (Ix, 2x and 4xRT, 
respectively) 



Power spectrum onlyEI Power spectrum only 



Usability Numerical simulations Numerical or semi- Numeiical or semi- Analytical model- 

numerical simulations numerical simulations ing (no realization), 

numerical or semi- 
numerical simulations 



Well defined 




Yes 


Inaccurate assumptions 
small scales 


on 


Unable to use on a grid too 
fine (see 516.5.31 


Yes 


Error0 in 
ID 


at 

k ^ 2/i/Mpc 


< 2% 


benchmark 


<4% 


< 2% 


0% 


< 1% 


0% 


< 10% 


power 
spectrum 


at 

2<k<7h/Mpc 


< 20% 


benchmark 


< 40% 


< 20% 


< 1% 


< 14% 


< 5% 


< 10% 


SUs 


Preliminarjifl 








350 


358 


375 


350 


358 


375 


(=cores 


Processing 


2048 


2127 


0.1 


0.7 


8.5 


52 


887 


5.3 


x hours) 


Total 


2048 


2127 


350 


359 


384 


402 


1245 


380 



In principle, a brightness temperature data cube in redshift space can be constructed by taking the inverse Fourier transform of the k-space biightness 
temperature evaluated using equations l |57t and l |59t for the DEMRF scheme and quasi-linear /ij^-decomposition scheme, respectively. However, ahasing 
effects from multiple forward and backward FFTs can introduce errors. It is beyond the scope of this paper to test these effects. 

All errors here are with respect to the results from the PPM-RRM (4xRT) scheme, which is the most accurate. 

Preliminary SUs for the MM-RRM scheme, quasi-linear /^k-decomposition scheme, and DEMRF scheme refers to the SUs used to smooth particle density 
and velocity data onto a regular, real-space, giid. 



6.7 Computational Accuracy and Efficiency 

So far we have discussed four viable schemes to compute 21cm 
brightness temperatures in redshift space: the PPM-RRM scheme 
(§ 16.2.1b . the MM-RRM scheme (§ I6.2.2t . the quasi-linear ^k- 
decomposition scheme (§|533J, and the DEMRF scheme (§|633]l. 
To facilitate the usage of these schemes, we compare their usability, 
accuracy and efficiency in Table[T] 

With N-body particle data (numerical simulation), the PPM- 
RRM scheme has no ambiguity in finding new particle locations in 
redshift space. When particle data is re-smoothed onto a redshift- 
space grid four times finer than RT grid resolution, PPM-RRM 
(4xRT) can accurately compute the power spectrum down to the 
RT resolution scale. However, the scheme is very computationally 
expensive and difficult to code, so we recommend to use it only as 
a development tool and for benchmarking, not for production work. 

The MM-RRM (4xRT) scheme is the perfect tool for produc- 
tion work. It requires only 1/6 of total computational effort of the 
PPM-RRM (4xRT) scheme (including preliminary calculations), 
and the results are just as accurate. Using the fine (4xRT), instead 
of coarse (RT) grid does not really add to the total computational 



effort. Note also that it can be directly used for semi-numerical sim- 
ulations that evolve density on grids and do not use particles. 

The DEMRF scheme is a nicely posed scheme since it is just 
a mathematical integration. However, in practise if we wish to sub- 
stitute the cell-wise average ^exp ^—i (^jij§^^^ '^ll'^ll ^ 

with exp i ^ °°.°) ^ '^H (^11 (''))coii] '^^''^8 cell-wise velocity, 
this scheme loses accuracy at the cell size ~ 114/1024 ~ 
0.11 Mpc/h. Moreover, the DEMRF (2xRT) scheme is neither the 
most accurate nor the most efficient, so it is not to be recommended 
for production work. However, it is useful for validating the results 
from the PPM-RRM and MM-RRM schemes. 

In the case of no realization, one can employ the quasi-linear 
^k-decomposition scheme which yields the redshift-space power 
spectrum with moderate accuracy and the least computational ef- 
fort. It also has as a useful feature that it can proceed with only 
real-space statistics as input, making it an ideal tool for pure ana- 
lytical modeling. 

The upshot is that we recommend the MM-RRM (4xRT) 
scheme for practical usage, and the PPM-RRM (4xRT) for devel- 
opment. 
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Figure 8. Test of the MM-RRM scheme: 21cm redshift-space ID power 
spectrum at z = 9.457 (50% ionized). We experiment with grids of the 
RT grid resolution (short-dashed, green), 2 times (long-dashed, blue), and 
4 times higher resolution (dot-short-dashed, cyan). The benchmark is the 
result from the PPM-RRM (4xRT) scheme (solid, black). The power spec- 
tra of the MM-RRM (4xRT) scheme and the PPM-RRM (4xRT) scheme 
are indistinguishable on all scales shown. The fractional en'ors of the MM- 
RRM results with respect to the PPM-RRM result are shown in the inset. 



7 HOW MUCH DO RARE OPTICALLY-THICK CELLS 
AFFECT THE ACCURACY OF POWER SPECTRUM IN 
THE OPTICALLY-THIN APPROXIMATION? 

In the optically thin limit, computation of 21cm redshift-space 
brightness temperature can be simplified by taking advantage of the 
proportionality of brightness temperature and neutral atom density, 
both in redshift-space. However, we have shown in Figure |3] that 
there is a nonzero, albeit small, chance to find large 21cm optical 
depth in the IGM. So in principle, the observed 21cm power spec- 
trum in the redshift-space that takes the finite optical depth into ac- 
count can be different from the result in the optically-thin approxi- 
mation. The difference depends on the population of optically-thick 
cells. In this section, we revisit the accuracy of the optically-thin 
approximation with regard to the 21cm power spectrum. 

We use the DEMRF method to calculate two power spectra in 
redshift-space: one with finite optical depth in equation i56\ , and 
one in the optically-thin approximation in equation ( I57t . We have 
demonstrated in i; 16. 5. 3 1 that, in the optically thin limit, the power 
spectrum using the DEMRF scheme agrees with the PPM-RRM re- 
sult in the comfort zone (k ^ fcjv/4). Here, we smooth the density, 
velocity and velocity gradient fields on the fine (512'') grid with 
two times better resolution than RT grid (256^), and focus on the 
region k ^ fc^^^V4 = 3.5 Mpc/h. We use SPH-like smoothing 
of our N-body particle data to compute the velocity gradient on the 
grid. Details of this technique are discussed in Appendix IA] 

The 21cm optical depth depends on the spin temperature 
which, however, is beyond the scope of our cosmological radiative 
transfer simulations. For the purpose of demonstration, we assume 
a = Ts /TcMB (zcos ) is a spatial constant, and investigate the cases 



Q — 100, 10, and 0.1 (the a = 1 case has no 21cm radiation con- 
trast to CMB). In the optically-thin approximation, the power spec- 
trum with finite (but constant) spin temperature is just the power 
spectrum with high spin temperature (Ts 2> Tcmb), i.e. the result 
in ^ 16331 scaled by the factor {1- ^f. 

In Figure |9] (left panel), we find that the power spectra in the 
optically-thin approximation are so highly accurate, in the a = 100 
and 10 cases, that the two curves (finite optical depth vs. optically- 
thin approximation) are almost indistinguishable. However, Fig- 
ure |9] (right panel) shows that, in the low Ts case (a = 0.1), the 
optically-thin approximation can result in an error of ~ 10% in 
the power spectrum on large scales, and > 30% on small scales. 
The large-scale error is due to the offset in the global mean sig- 
nal, because the optically-thin approximation overestimates the 
brightness temperature (i.e. STt oc Tj,) in the optically-thick cells, 
which should otherwise be suppressed in the exact form STt oc 
[1 — exp (— T^)] when optical depth is large. This decreases the 
small-scale power spectrum, too, because the 21cm brightness tem- 
perature in these overdense regions (where Ti, > 1) fails to encode 
the complete statistical information of density and ionization fluc- 
tuations. 

Is the optically-thin approximation accurate with regard to 
21cm power spectrum? The answer depends on the spin temper- 
ature, because 21cm optical depth is inversely proportional to Ts- 
As Figure |3] shows, the low Ts case has much higher chance to 
find optically-thick cells than the high Ts case, i.e. roughly an or- 
der of magnitude smaller in Ts, an order of magnitude larger in the 
probability of > 1. This is consistent with our results that the 
optically-thick cells are too rare to virtually affect the power spec- 
trum when Ts/TcMB ^ 10, but they are non-negligible when Ts 
is lower than Tcmb . 

The upshot is that the power spectrum in the redshift-space 
calculated in the optically-thin approximation, e.g. using the PPM- 
RRM or MM-RRM scheme, is accurate with respect to the result 
that takes finite optical depth into account, only when Ts is high 
{Ts/TcMB ^ 10). The low Ts case merits further careful investi- 
gation which we defer to future work. 



8 HOW ACCURATE IS LINEAR THEORY? 

The l inear theory formula f or 21cm redshift-space power spec- 
trum dSarkana & Loeb 2005) ha s been widely employed in the lit- 
erature (e.g. ISantos & Coorav. ,2006t IZahn et al.ll2007l: iMao et al.l 
I2OO8I : lAdshead et al.l 201 1 ), but it is derived under two assump- 
tions that may both break down. First, the ionization fluctuations 
are assumed to be linear. This is only valid on scales much larger 
than the size o f the H II region whi ch can be rather large (~ 10 
Mpc, see, e.g. iFriedrich et al]|201ll) . Second, the matter density 
and velocity fluctuations are assumed to be linear, i.e. the veloc- 
ity is dictated by the density through the linear relation, ^[[^(k) = 

' H(Zco.) ' _ 
fe 



1+z ) ^phC^)^- '^^i^ relation is also inaccurate on small 
scales. Is linear theory spoiled by the breakdown of these approx- 
imations? For simplicity, we restrict our discussion in this section 
to the simple case Ts S> Tcmb ■ 

We compute the 21cm redshift-space ID power spectrum in 
the linear theory by angle-averaging equation ( I60t with moments 
in equations - {69}, tMcOuinn et al..,2006, ; ,Zahn et al.„2007l : 
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Figure 9. Power spectra of 21cm brightness temperature in redshift-space calculated in the optically-thin approximation (dotted, blue), and the results that 
take finite optical depth into account (long-dashed, red), both using the DEMRF scheme on a grid (512^) two times finer than the RT grid. (Left) When Ts 
is high, i.e. a = Ts/Tqmb = 100 (thin lines) and 10 (thick lines). In each set, two curves (finite optical depth vs. optically-thin approximation) overlap 
almost exactly. (Right) When Ts is low, i.e. a = 0.1. All results use the RT simulation data at 50% ionized epoch (2 = 9.457). The fractional errors of the 
optically-thin approximation are with respect to the results with finite optical depth. The vertical lines at A: = fc^^^'/4 = 3.5 h/Mpc delimit the comfort zone 
for the D£MSF results computed on the 512"^ grid. 



iLidz et al.ll2007h 



\k) = ST, 



PL 



3' ''-hi'^Vh 



+—P^ s (k) 



(72) 



We compare it with the angle-averaged fully nonlinear power spec- 
trum in redshift-space, computed using the PPM-RRM (4xRT) 
scheme. In Figure (TO) we find that the linear theory power spec- 
trum departs from the fully nonlinear result with < 30% error in 
the intermediate range of fc ~ 0.1 — 1 ft/Mpc, at the 50% ionized 
epoch. It crosses the nonlinear result at A; ~ 1 /i/Mpc, and devi- 
ates more from the latter at smaller scale s. This is in qua litative 
agreement PI with a similar comparison in lLidz et al.l (1200 7) (their 
Fig. 10, but they did not provide any detail of how they computed 
the nonlinear power sp ectrum in redshift-space). 

iLidz et all ( |2007 |) pointed out that such a large error in linear 
theory may result from the neglect of higher order auto- and cross- 
correlations involving density and ionization fluctuations, i.e. the 
breakdown of the first assumption we mentioned above, but they 
did not provide a solution that incorporates all of relevant higher 
order terms in redshift-space power spectrum. Here we propose in 



^'^ Th e deviation increases monotonically at k > Ih /Mpc in iLidz et al] 
while there seems to be a tum-around at large k in our Figure 1101 
This turnaround is not real because it is the resolution effect. We are forced 
to compute the linear theory power spectrum on the RT grid (the grid for 
ionization fraction fields). The aliasing effect suppresses the linear theory 
result at k > k'-^^^^ /i = 1.76 h/Mpc in our simulation, while Lidz et al. 

result is free of aliasing effect at fc < 10 /i/Mpc by adopting an RT 
grid of much higher resolution. 



§[533]the quasi-linear ^^^-decomposition scheme as such a solu- 
tion that not only incoiporates these higher order corrections, but 
can decompose 2Icm redshift-space power spectrum in polynomi- 
als of /ik, just as the linear theory does. How accurate is this new 
scheme? Figure [TO] also shows that the angle-average power spec- 
trum of the quasi-linear /ik -decomposition scheme (calculated us- 
ing eq. [5j agrees with the fully nonlinear result to ~ 10% accu- 
racy at ~ 0.3 — 2 /i/Mpc, but with increasing errors at smaller 
scales. We will defer the detailed investigation of the errors in this 
scheme associated with the neglect of additional nonlinearity to 
IShaDiroetallfeOllI) . 

The large errors in the linear theory for redshift-space distor- 
tion suggest that it is a simple, but by no means accurate, tool to 
predict 21cm power spectrum. One should either employ the quasi- 
linear iiu,-decomposition scheme for improved (but not perfect) ac- 
curacy, or follow the numerical schemes we proposed above (PPM- 
RRM and MM-RRM) to obtain fully nonhnear results. 



9 HOW ACCURATE IS THE "Vu-LIMITED" 
PRESCRIPTION? 

9.1 The "VD-limited" prescription vs. the avoidance of the 
divergence problem in observer space 

ISantos et al. treated the effects of peculiar velocity on the 

21cm brightness temperature by evaluating an equation equivalent 
to our equation ( I36t at each point in a real-space grid at a given 
time. They found that the 21cm brightness temperature diverges 
in some overdense regions where Sg^^ — >■ —1. As such, the power 
spectrum computed from the Fourier transform of this 21cm bright- 
ness temperature evaluated in real-space diverges, too. They deal 
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Figure 10. Test of the linear theory: 21cm redshift-space ID power spec- 
tr um at 2 = 9-457 (50% ionized), calculated using the linear theory 
of ISarkana & Loebl feOOSi) (dot-short-dashed, blue), the quasi-lineai' /x^- 
decomposition scheme (long-dashed, red), and the PPM-RRM (4xRT) 
scheme (solid, black), respectively. 



with this divergence problem by replacing the actual value calcu- 
lated for Sg^^ from their real-space grid data whenever it is close 
to -1, by a fixed minimum value larger than —1 (e.g., —0.7 in their 
paper), so as to cap the divergence and obtain finite results for both 
brightness temperature and its power spectrum. This approach was 
also a dopted by the llcmFAST code l lMesinger. Furlanetto & CenI 
1201 ih (with the cap -0.5). 

Before analyzing the accuracy of the V«-limited prescrip- 
tion, we would like to explain why the divergence encountered 
for &g^y —1 is a mathematical, but not a physical one. As we 
shall show, the appearance of the divergence is avoided naturally 
for physical observables in observer redshift-space. 

The first part of this explanation was already considered in 
§ 15.1.51 Equation ( 136b was derived under the assumption of low 
optical depth. However, the locations at which Sg^^ approaches - 
1 are not optically thin. The 21cm brightness temperature must be 
evaluated using equation ( I47l >. instead, at these locations, to take 
finite optical depth into account. When this is done, the brightness 
temperature does not diverge for Sg^^ — >■ — 1. 

However, even in the optically-thin approximation, it is un- 
necessary to apply a cap to the velocity gradient in order to pre- 
vent divergence in the physical observables, as long as we ac- 
count properly for redshift-space distortion. The approach in which 
equation ( 136b is applied to real-space grid data does not fully ac- 
count for the remapping of real- to redshift-space locations of 
21cm sources. While this remapping cannot remove the diver- 
gence of 21cm brightness temperature at those locations at which 
5g^v — >■ — 1, the power spectrum in redshift-space is guaranteed 
to be finite. The reason is simply that real-space regions for which 
Sg^v — >■ — 1 become infinitesimally small in redshift-space, since 
dp s = dPr j 1 + (5g^^ (r) | . The Fourier transform of brightness tem- 
perature in redshift-space, 5T§{\i) = / d^s e^'*' "' 5T§(s), is a fi- 



nite integration, and so is the power spectrum computed from it, 
because the divergent factor in 5T§{s,) — 5Tl{v) oc 1/ |l + Sg^^\ 
is exactly cancelled by its inverse in the volume element dPs . In 
addition, in the optically-thin approximation, while, strictly speak- 
ing, the 21cm brightness temperature still diverges at those loca- 
tions at which 5g^,^ approaches -1, it, too, becomes finite when 
smoothed over finite band- and beam-width in observer redshift 
space (see also ij |5.2.3b . This, of course, makes perfect sense since 
the 21cm emitted photons produced by a finite region of space, in 
the optically-thin limit, are proportional to the number of neutral 
hydrogen atoms in that region, which is always finite and is con- 
served by the mapping from real- to redshift-space. Since what ob- 
servers actually measure are this pixelized brightness temperature 
and the power spectrum, full account of redshift-space distortion 
gives a physical result for these observables without resorting to 
artificial caps. 

9.2 Evaluating the accuracy of the "Vu-limited" 
prescription 

Although based on a conceptual artifact (truncation of an un- 
physical divergence) and providing an incomplete fix (calculat- 
ing the power spectrum in real- instead of redshift-space), the 
"Vw-limited" prescription may still provide a practic al solution to 
the pr oblem of the diverging brightness temperature. ISantos et al.l 
( |2010 b argued that although an ad-hoc solution, imposing this limit 
only affects a very small number of cells, and thus has no influence 
on global statistics such as the power spectrum. Since our methods 
avoid the divergence problem, we are now able to test this asser- 
tion. Furthermore, to be a practical solution to the problem, the 
results should not depend too much on the choice for the cap on 
the velocity gradient. Here we test these two issues by comparing 
the results of the V?;-limited prescription to those from our PPM- 
RRM (4xRT) scheme. For simplicity, we restrict our discussion in 
this section to the simple case Ts 3> Tcmb • 

In order to find gridded values for the velocity gradient, we 
use SPH-like smoothing of our N-body particle data to compute 
the velocity gradient. Details of this technique are discussed in Ap- 
pendix [A] We implement the Vu-limited prescription by replacing 
the actual value of {Sg^^)^^^-^ by the cap value of —A whenever 
{5g^v)^^ii < —A, for a range of cap values A =0.1, 0.3, 0.5, 0.7, 
and 0.9, evaluating and Fourier transforming the brightness temper- 
ature in real space. We then average the power spectra over three 
LOS directions. In this section, we assume the limit of high spin 
temperature, Ts » Tcmb- 

To most clearly show the effects of the Vn-limited prescrip- 
tion we first take our volume to be fully neutral, by setting xm = 1 
everywhere. Figure [TT] (top left panel) shows the power spectra 
from the Vu-limited prescription for five different values of the 
upper limit A as well as the power spectrum calculated with the 
PPM-RRM scheme. Here we use the smoothed velocity gradient 
field on the RT grid resolution. We find that different values for A 
yield rather different power spectra even on large scales, and none 
of the previously proposed values of caps (A — 0.5 or 0.7) is con- 
sistent with the PPM-RRM result. 

This is of course the most extreme case since in a fully neu- 
tral medium all locations with Sg^^ — >■ —1 contribute. Since 
these regions are preferably located in high density areas, which 
typically reionize earlier, one can expect that the effect is much 
less severe when considering a neutral fraction distribution ob- 
tained from a reionization calculation. Figure [TT] (top right panel) 
shows this indeed to be case. On large scales, different limits 
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Figure 11. Power spectra of 21cm brightness temperature at 50% ionized epocli (z = 9.457). We sliow the results of the Vti-limited prescription with the 
cap {^drv)i,i,ii is ~^ at A =0.1 (dotted, red), 0.3 (short-dashed, brown), 0.5 (long-dashed, blue), 0.7 (dot-short-dashed, cyan), and 0.9 (dot-long-dashed, 
magenta), compared to the results of the PPM-RRM (4xRT) scheme (solid, black). The fractional errors of the Vu-limited prescription are with respect to the 
PPM-RRM (4xRT) result, plotted in the inset. Upper panels: using velocity gradient field smoothed on the RT grid resolution; lower panels: using velocity 
gradient field smoothed on the fine grid (4 times finer than RT grid resolution). Left panels: assuming a fully neutral Universe (i'hi = 1); right panels: using 
the actual reionization fluctuations from the simulation. 



in the Vu-limited prescription yield converging power spectra 
whicli, liowever, have an offset of ~ 20% from the power spec- 
trum of the PPM-RRM scheme. This offset is due to the en- 
hancement in the mean brightness temperature averaged in real 
space, i.e. although the distribution of (5s^i,(r) has zero mean, 
the distribution of (1 + (5pHi(r)) / |1 + Sa^,v{'r)\ does not have the 
(volume-weighted) mean of unity [3 due to the nonlinear func- 



This can be compared to the mean in redshift space, where the averaging 



tion 1/ |1 + 5g^y\. Although converged on large scales, on small 
scales the results of the Vu-limited prescription still depend on the 
cap value chosen and can have inaccuracy as large as < 40% for 
A = 0.5, or < 50% for A = 0.7, both at fc < 2/i/Mpc, and more 



is over redshift-space volume elements d"^ s = (fir \1 + Sg equiva- 
lent to averaging in real-space weighted by 1 1 + (5g „ (r) | , and therefore the 
distribution of (1 + <5pjjj (r)) / 1 1 -f ^s^i, (r)] has 1 as the mean in redshift 
space. 
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divergent for larger caps (as exemplified by A = 0.9). However, 
these inaccuracies are substantially smaller than the ones found for 
the fully neutral case. 

These inaccuracies, of course, depend on the grid resolution 
of the velocity gradient field. We redo the above analyses, using a 
fine grid (four times finer than the RT grid), as shown in Figure [TT] 
(bottom panels). We find that the errors in the results of the Vv- 
limited prescription are significantly amplified. This is because the 
velocity and its gradient become more nonlinear on smaller scales. 
Hence, a larger population of cells are "clipped" on fine grids than 
on coarse grids. 

The reasons why the Vu-limited prescription does not work 
well are twofold. First, this prescription deals with data cubes in 
real-space coordinates. Consequently, their Fourier transform and 
power spectra are real-space quantities, unlike in the redshift-space 
as this prescription claimed to do. Second, even though the Vv- 
limited prescription was invented to circumvent the unphysical di- 
vergence in 21cm brightness temperature in real-space regions that 
are optically thick to 21cm radiation, the cells that are affected by 
imposing a cap on velocity gradient are actually much more nu- 
merous than the optically-thick cells; e.g., at 50% ionized epoch 
in our simulation, we find that only a fraction of ~ 10^^ amongst 
all cells are optically thick in the best case (Ts/Tcmb = 100), or 
~ 0.01% in the worst case (T^/Tcmb = 0.1), (see ^ |5X6] l, but 
the Vw-limited prescription affects all those cells for which \5g^v\ 
exceeds the cap, which can be a much larger fraction of the cells 
than that of the optically thick cells. A fraction ~ 1% of the cells 
have Sg^v ^ —0.7, while ~ 3% have Sg^y ^ —0.5. (These spe- 
cific fractions can depend on the grid resolution of the simulation. 
The numbers here are counted using velocity gradient field on RT 
grid resolution.) In other words, the Vf-limited prescription affects 
a lot more cells than necessary. 

Even though the V?;-limited prescription cannot provide the 
most accurate treatment of the effect of peculiar velocity, it still 
serves as a simple and useful tool to estimate the 21cm power spec- 
trum, if a smaller cap is chosen than those previously proposed. We 
optimize this prescription by comparing its results using A — 0.1 
and 0.3 with our PPM-RRM result, and find that with the actual 
reionization fluctuations, the Vf-limited prescription with A = 0.1 
approximates the PPM-RRM result with the least errors < 20%, 
while, if we assume a fully neutral universe, A = 0.3 is the optimal 
cap, with errors < 10%. The optimal value of the cap depends on 
the grid resolution, and perhaps on the redshift and the ionization 
fraction, as well. 



10 CONCLUSIONS 

• We have demonstrated that the neglect of peculiar velocity in- 
troduces a substantial error in 21cm brightness temperature spec- 
tra from the FOR and noticeable anisotropy in the 21cm power 
spectrum. We did this in three different ways: first, we compared 
the 3D power spectra computed uncorrected for peculiar velocity 
(UPV scheme), from the quasi-linear /Xk-decomposition scheme, 
and from a particle-based numerical scheme (PPM-RRM); second, 
we compared the 21cm brightness temperature spectra computed 
from the UPV scheme and the PPM-RRM scheme, along 5 differ- 
ent sightlines; lastly, we compared the angle-averaged 21cm power 
spectra computed from the quasi-linear ^k-decomposition scheme 
and the UPV scheme, respectively. The non-trivial difference be- 
tween results with and without peculiar velocity correction moti- 
vates our thorough investigation of the effect of peculiar velocity 



on 21cm signal as well as a more careful treatment of this effect on 
reionization simulation data than previously made. 

• We clarify that peculiar velocity distorts the mapping of 21cm 
brightness temperature not only by shifting the apparent location in 
redshift-space, but also by modifying the brightness temperature it- 
self in real-space. We show that the combined effect, which we call 
"21cm redshift-space distortion", establishes, in the limit of low 
optical depth and high spin temperature, the exact proportionality 
between observed 21cm brightness temperature and the neutral hy- 
drogen density as measured in redshift-space. This proportionality 
makes it possible to infer the three-dimensional distribution of neu- 
tral hydrogen density using 21cm brightness temperature measure- 
ments. 

• We show that this proportionality between 21cm observed 
brightness temperature and the redshift-space neutral hydrogen 
density, however, can break down when T2icm 1 and/or Ts < 
TcMB- For the first case, we check the optically thin approxima- 
tion, and demonstrate that this widely-assumed approximation is 
mostly valid in the IGM, but it can be invalid in some cases, e.g. 
in virialized halos where the peculiar velocity gradient can be large 
enough to cancel the Hubble flow. For the second case, we show 
that the proportionality mentioned above is spoiled by the spatially- 
varying r,,-dependent factor 1 - Tcmb/TJ' (r). This TJ'°^(r) 
includes a correction to Ta of the order v/c due to peculiar velocity. 

• The unphysical divergence in 21cm brightness temperature re- 
sults from the neglect of finite optical depth, which eliminates the 
divergence. We show that, in the optically thick limit, the optical 
depth can depend upon higher order spatial derivatives of peculiar 
velocity than dv^^ /dr^\ . 

• We derive the fully nonlinear Fourier transform of 21cm 
brightness temperature fluctuations, with finite optical depth, as 
measured in redshift-space, in terms of the density, velocity and 
its gradient, ionization fraction, and spin temperature fields in real 
space, following the combined effect of 21cm redshift-space dis- 
tortion. We further simplify it in the optically-thin approximation. 
We further show that, when redshift-space distortion is properly 
accounted for, however, the observed power spectrum in redshift- 
space remains finite even in the optically-thin approximation. 

• We investigate the effect of finite 21cm optical depth. The 
21cm power spectrum in redshift-space calculated in the optically- 
thin approximation is accurate with respect to the results which take 
finite optical depth into account, only when spin temperature is high 
relative to the CMB temperature (Ts/Tcmb > 10). 

• We clarify that it is the bulk velocity of the gas but not the 
thermal velocity that is responsible for the velocity correction to 
the optical depth and 21cm brightness temperature. This is done by 
showing that the latter constitutes only a negligible contribution to 
the correction, compared to the former, when r2icm < 1. 

• To make a careful treatment of the peculiar velocity effect on 
21cm brightness temperature when using reionization simulation 
data, we propose and test two numerical schemes that compute 
the 21cm brightness temperature as measured in a redshift-space 
grid from real-space simulation data, in the limit of high spin tem- 
perature. Both schemes take advantage of the mapping from real- 
to redshift-space, one particle-based (PPM-RRM), and one grid- 
based (MM-RRM). We show that the MM-RRM scheme can be 
optimized to achieve the same high accuracy in the angle-averaged 
power spectrum as the PPM-RRM scheme, while being much more 
computationally efficient than the latter. If the RT grid resolution 
(on a mesh with Nyquist wavenumber fcjv,RT, the mesh on which 
the ionization fluctuation field is determined) is coarser than the 
resolution of the density and peculiar velocity fields, we optimize 
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the grid-based MM-RRM resolution by including all modes with 
k ^ fcjv^RT in a grid with Nyquist wavenumber 4fcjv,RT which 
uses the finer-resolution density and velocity data, together with 
the coarser-resolution ionized fractions. This reduces the aliasing 
errors which would otherwise spoil the results for k > fcjv,RT /4 
if all data were coarsened to the RT-grid resolution. We show that 
this optimized MM-RRM scheme can compute the angle-averaged 
21cm power spectrum within < 1% error with respect to the PPM- 
RRM(4xRT) results, for all modes fc < fcjv.RT- 

• We examine the linear theory formula widely em- 
ployed to compute th e 21cm redshift-space power spectrum 
( lBarkana&Loebll2005h . and find large inaccuracy (~ 30%) at 
the intermediate range k ~ 0.1 — 1 /i/Mpc at the 50% ionized 
epoch, in the high spin temperature regime. This suggests that 
linear theory cannot work as an accurate tool to predict the 21cm 
power spectrum in redshift-space. 

• We develop the "quasi-linear /Xk-decomposition scheme" 
which can decompose 21cm power spectrum in polynomials of /^k, 
just as the linear theory does, but it incorporates relevant higher or- 
der correlations of ionization and density fluctuations. We find that 
the fully nonlinear 21cm ID power spectrum deviates from the pre- 
diction of quasi-linear /ik -decomposition scheme by roughly 10% 
at the 50% ionized epoch (see i; 16.5.21 ). The nonlinearity may intro- 
duce larger deviations when the 3D power spectrum is decomposed 
to extract only the P^4, (k) for cosmology. It is important to under- 
stand the nature of this nonlinear effect, and estimate its impact on 
21cm cosmol ogy. We will addres s these issues in the second paper 
of this series dShapiro et al.l201l[) . 

• Our careful treatment of brightness temperature fluctuations in 
redshift space avoids the divergences that appear in the real-space 
evaluation when peculiar-velocity gradients are large. Such large 
gradients are a natural result of nonlinear structure formation on 
small scales. We find that previous attempts to escape these diver- 
gences by numerically "clipping" the velocity gradients whenever 
they exceed some threshold (refeiTed to here as the "Vw-limited 
prescription") introduce a non-negligible inaccuracy in the 21cm 
power spectra on all scales, including scales much larger than that 
of the nonlinearity. We show that the errors associated with this pre- 
scription, however, can be reduced if the value of the cap is properly 
chosen (e.g. A ~ 0.1 yields an error ~ 15% at fc ~ 0.1 fe/Mpc), 
but this error grows with increasing spatial resolution of the grid, 
and may depend on redshifts and ionization fraction, too. 

• The upshot is that we provide an integrated understanding of 
how peculiar velocity affects 21cm tomography, and also an ac- 
curate and efficient numerical algorithm (MM-RRM) for practical 
numerical application. 
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(j = 1, . . . , Np). We define a particle's kemel hi to be the distance 
between the particle i and its 32""* nearest neighbor particle. We 
take the "scatter" a pproach to smooth particle data (see Fig. 2 of 
IShaDiroetal.lll996l for an illustration of the scatter vs. gather ap- 
proaches), i.e., a field point at r is influenced by a particle i if this 
particle's own influence zone covers this field point (e.g., in the case 
of isotropic kemel, |r — r^] ^ hi). 

We employ the triangular kemel function with adaptive kemel 
size h, W{v; h) = fh(x)fhiy)fh{z), centered at the particle lo- 
cation to smooth its data. The function fh{x) is triangular-shaped 
with width 2h, i.e., 

f + i , O^x^h 

.fh{x) = i ^ + i , ^h^x<0 (Al) 
^ , otherwise 

Smoothed Fields at a Point 

The smoothed mass density and momentum density fields are de- 
fined, respectively, by 



p(r) = y^^miWjr ~ ri;hi) , 

i 



(A2) 
(A3) 



To preserve momentum, the continuous velocity field is defined by 

v(r) =P(r)/p(r). (A4) 

We identify the bulk-flow velocity of the IGM at a particle's posi- 
tion Ti to be the smooth field v(ri) evaluated at r^. 



Smoothed Fields of a Cell 

To smooth particle data onto a regular grid, we use the following 
approach to compute the cell-wise mass density, 



''cell Jccll 



TT—y^mi W{r ~ ri;hi)d^r , 

VccU — ' ' ■■ 



(A5) 



where the integral J^^jj W{r — r^; hi)d?r can be evaluated analyti- 
cally, and is only a function of hi and the relative location between 
the particle i and the cell boundaries. Similarly, the cell-wise mo- 
mentum density is 

(^)ceii = / n^)d'r 

Vccll Jccll 

= -i-Vm^v, / W{v~r,-h.)d\. (A6) 

Vccll ^ Jccll 

The cell-wise velocity is defined in a momentum-preserving way, 

(V>ccll = (^)ccll/(p)ccll • (A7) 



APPENDIX A: SPH-LIKE SMOOTHING WITH 
ADAPTIVE KERNEL 

In this section, we briefly describe the SPH-like technique to 
smooth N-body pa rticle data onto a grid. We refer readers to 
IShapiroetal1 ( ll996h for a comprehensive discussion of smooth par- 
ticle hydrodynamics with an adaptive kemel. 

Assume that the continuous density and velocity fields are rep- 
resented by A'^p particles with mass rrii, location r^, and velocity Vi 



dv |/dr | of a Cell 

We compute the cell-wise velocity gradient (dun /dr|| 
following way (assuming the LOS is along one of the principal 
axes of a cubical cell), 

dv\\ 



(^dv\\/dr 



Vccii Jccll rfj"!! 

at[( 



(r) 

\ II / — plaiicj 



plane 



(A8) 
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Figure Al. Caitoon of computing (^dv^^/dr^^ ^^^^^ for a cubical cell of size 
AL on each side. In this cartoon, we assume the LOS along the xs-axis, 
then the "+ plane" ("— plane") is the X1-X2 plane with X3 = AL {x = 0). 



where AL is the size of the cubical cell, "+ plane" ("- 
plane") is the cell wall perpendicular to the LOS with larger 
(smaller) location along the r||-axis, and ('"|| )_|, pij^^^, is the 
velocity mean on the "+" cell wall, i.e. (v«) , , — 

\ II / + plane 

j^^\yA /+ piano '^ll ('S'ii + Pl^-iis). Unfortunately, we cannot 
apply the same smoothing as in equation JA5I > to compute the ve- 
locity average, because the velocity defined in equation ( IA4b in- 
volves a summation in the denominator 

To circumvent this, we approximate the smoothed velocity av- 
eraging on a cell wall by the momentum-preserving velocity, i.e. 

('^Il)planc^<^ll)plane/(P) plane, (A9) 

where the rh.s. is the center-of-mass velocity of a thin layer on the 
cell wall. The surface mass density of the cell wall is 

(P>planc = 7T7\2 / P(^)'^^ 



(AL)2 



plane 



■E 



plane 



(AL)2 

where the integral can be evaluated analytically. 



; ftOrf 7-x , (AlO) 



Wir-r^-h,)d\_L 



plane 



;i_c + AL/2 



-AL/2 



fh,{xi - Xl^i)dxi 



X2,^+AL/2 



a:2,c-AL/2 



fh,{x2 - X2,i)dX2 



X fh, (a;3,pla 



■ X3,, 



(All) 



Here we take xi and X2 to be axes in the cell wall perpendicular 
to the LOS axis X3, 13, piano is the LOS coordinate of the cell wall, 
Xi^c and 0:2, c are the transverse coordinates of the center of the cell, 
and Xi are the three-dimensional coordinates of the particle i. 
Similarly, we use 



1 



plane 



(Ai)- ./plane 



P||(r)dVx 



■E 



plane 



W{r-ri;hi)d''r^. 



